OpenVDB  2.3.0
VolumeToMesh.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2013 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 
31 #ifndef OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
32 #define OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
33 
34 #include <openvdb/Platform.h> // for OPENVDB_HAS_CXX11
35 #include <openvdb/tree/ValueAccessor.h>
36 #include <openvdb/util/Util.h> // for COORD_OFFSETS
37 #include <openvdb/math/Operators.h> // for ISGradient
38 #include <openvdb/tools/Morphology.h> // for dilateVoxels()
39 #include <openvdb/tree/LeafManager.h>
40 
41 #include <boost/scoped_array.hpp>
42 #include <boost/scoped_ptr.hpp>
43 #include <tbb/blocked_range.h>
44 #include <tbb/parallel_for.h>
45 #include <tbb/parallel_reduce.h>
46 
47 #include <vector>
48 #include <memory> // for auto_ptr/unique_ptr
49 
50 
52 
53 
54 namespace openvdb {
56 namespace OPENVDB_VERSION_NAME {
57 namespace tools {
58 
59 
61 
62 
63 // Wrapper functions for the VolumeToMesh converter
64 
65 
72 template<typename GridType>
73 void
75  const GridType& grid,
76  std::vector<Vec3s>& points,
77  std::vector<Vec4I>& quads,
78  double isovalue = 0.0);
79 
80 
89 template<typename GridType>
90 void
92  const GridType& grid,
93  std::vector<Vec3s>& points,
94  std::vector<Vec3I>& triangles,
95  std::vector<Vec4I>& quads,
96  double isovalue = 0.0,
97  double adaptivity = 0.0);
98 
99 
101 
102 
105 
106 
109 {
110 public:
111 
112  inline PolygonPool();
113  inline PolygonPool(const size_t numQuads, const size_t numTriangles);
114 
115 
116  inline void copy(const PolygonPool& rhs);
117 
118  inline void resetQuads(size_t size);
119  inline void clearQuads();
120 
121  inline void resetTriangles(size_t size);
122  inline void clearTriangles();
123 
124 
125  // polygon accessor methods
126 
127  const size_t& numQuads() const { return mNumQuads; }
128 
129  openvdb::Vec4I& quad(size_t n) { return mQuads[n]; }
130  const openvdb::Vec4I& quad(size_t n) const { return mQuads[n]; }
131 
132 
133  const size_t& numTriangles() const { return mNumTriangles; }
134 
135  openvdb::Vec3I& triangle(size_t n) { return mTriangles[n]; }
136  const openvdb::Vec3I& triangle(size_t n) const { return mTriangles[n]; }
137 
138 
139  // polygon flags accessor methods
140 
141  char& quadFlags(size_t n) { return mQuadFlags[n]; }
142  const char& quadFlags(size_t n) const { return mQuadFlags[n]; }
143 
144  char& triangleFlags(size_t n) { return mTriangleFlags[n]; }
145  const char& triangleFlags(size_t n) const { return mTriangleFlags[n]; }
146 
147 
148  // reduce the polygon containers, n has to
149  // be smaller than the current container size.
150 
151  inline bool trimQuads(const size_t n, bool reallocate = false);
152  inline bool trimTrinagles(const size_t n, bool reallocate = false);
153 
154 private:
155  // disallow copy by assignment
156  void operator=(const PolygonPool&) {}
157 
158  size_t mNumQuads, mNumTriangles;
159  boost::scoped_array<openvdb::Vec4I> mQuads;
160  boost::scoped_array<openvdb::Vec3I> mTriangles;
161  boost::scoped_array<char> mQuadFlags, mTriangleFlags;
162 };
163 
164 
167 typedef boost::scoped_array<openvdb::Vec3s> PointList;
168 typedef boost::scoped_array<PolygonPool> PolygonPoolList;
170 
171 
173 
174 
177 {
178 public:
179 
182  VolumeToMesh(double isovalue = 0, double adaptivity = 0);
183 
184 
186 
187  // Mesh data accessors
188 
189  const size_t& pointListSize() const;
190  PointList& pointList();
191 
192  const size_t& polygonPoolListSize() const;
193  PolygonPoolList& polygonPoolList();
194  const PolygonPoolList& polygonPoolList() const;
195 
196  std::vector<unsigned char>& pointFlags();
197  const std::vector<unsigned char>& pointFlags() const;
198 
199 
201 
202 
205  template<typename GridT>
206  void operator()(const GridT&);
207 
208 
210 
211 
233  void setRefGrid(const GridBase::ConstPtr& grid, double secAdaptivity = 0);
234 
235 
239  void setSurfaceMask(const GridBase::ConstPtr& mask, bool invertMask = false);
240 
243  void setSpatialAdaptivity(const GridBase::ConstPtr& grid);
244 
245 
248  void setAdaptivityMask(const TreeBase::ConstPtr& tree);
249 
250 
254  void partition(unsigned partitions = 1, unsigned activePart = 0);
255 
256 private:
257 
258  PointList mPoints;
259  PolygonPoolList mPolygons;
260 
261  size_t mPointListSize, mSeamPointListSize, mPolygonPoolListSize;
262  double mIsovalue, mPrimAdaptivity, mSecAdaptivity;
263 
264  GridBase::ConstPtr mRefGrid, mSurfaceMaskGrid, mAdaptivityGrid;
265  TreeBase::ConstPtr mAdaptivityMaskTree;
266 
267  TreeBase::Ptr mRefSignTree, mRefIdxTree;
268 
269  bool mInvertSurfaceMask;
270  unsigned mPartitions, mActivePart;
271 
272  boost::scoped_array<uint32_t> mQuantizedSeamPoints;
273 
274  std::vector<unsigned char> mPointFlags;
275 };
276 
277 
279 
280 
288  const std::vector<Vec3d>& points,
289  const std::vector<Vec3d>& normals)
290 {
291  typedef math::Mat3d Mat3d;
292 
293  Vec3d avgPos(0.0);
294 
295  if (points.empty()) return avgPos;
296 
297  for (size_t n = 0, N = points.size(); n < N; ++n) {
298  avgPos += points[n];
299  }
300 
301  avgPos /= double(points.size());
302 
303  // Unique components of the 3x3 A^TA matrix, where A is
304  // the matrix of normals.
305  double m00=0,m01=0,m02=0,
306  m11=0,m12=0,
307  m22=0;
308 
309  // The rhs vector, A^Tb, where b = n dot p
310  Vec3d rhs(0.0);
311 
312  for (size_t n = 0, N = points.size(); n < N; ++n) {
313 
314  const Vec3d& n_ref = normals[n];
315 
316  // A^TA
317  m00 += n_ref[0] * n_ref[0]; // diagonal
318  m11 += n_ref[1] * n_ref[1];
319  m22 += n_ref[2] * n_ref[2];
320 
321  m01 += n_ref[0] * n_ref[1]; // Upper-tri
322  m02 += n_ref[0] * n_ref[2];
323  m12 += n_ref[1] * n_ref[2];
324 
325  // A^Tb (centered around the origin)
326  rhs += n_ref * n_ref.dot(points[n] - avgPos);
327  }
328 
329  Mat3d A(m00,m01,m02,
330  m01,m11,m12,
331  m02,m12,m22);
332 
333  /*
334  // Inverse
335  const double det = A.det();
336  if (det > 0.01) {
337  Mat3d A_inv = A.adjoint();
338  A_inv *= (1.0 / det);
339 
340  return avgPos + A_inv * rhs;
341  }
342  */
343 
344  // Compute the pseudo inverse
345 
346  math::Mat3d eigenVectors;
347  Vec3d eigenValues;
348 
349  diagonalizeSymmetricMatrix(A, eigenVectors, eigenValues, 300);
350 
351  Mat3d D = Mat3d::identity();
352 
353 
354  double tolerance = std::max(std::abs(eigenValues[0]), std::abs(eigenValues[1]));
355  tolerance = std::max(tolerance, std::abs(eigenValues[2]));
356  tolerance *= 0.01;
357 
358  int clamped = 0;
359  for (int i = 0; i < 3; ++i ) {
360  if (std::abs(eigenValues[i]) < tolerance) {
361  D[i][i] = 0.0;
362  ++clamped;
363  } else {
364  D[i][i] = 1.0 / eigenValues[i];
365  }
366  }
367 
368  // Assemble the pseudo inverse and calc. the intersection point
369  if (clamped < 3) {
370  Mat3d pseudoInv = eigenVectors * D * eigenVectors.transpose();
371  return avgPos + pseudoInv * rhs;
372  }
373 
374  return avgPos;
375 }
376 
377 
379 
380 
381 // Internal utility methods
382 namespace internal {
383 
384 template<typename T>
385 struct UniquePtr
386 {
387 #ifdef OPENVDB_HAS_CXX11
388  typedef std::unique_ptr<T> type;
389 #else
390  typedef std::auto_ptr<T> type;
391 #endif
392 };
393 
394 
396 enum { SIGNS = 0xFF, EDGES = 0xE00, INSIDE = 0x100,
397  XEDGE = 0x200, YEDGE = 0x400, ZEDGE = 0x800, SEAM = 0x1000};
398 
399 
401 const bool sAdaptable[256] = {
402  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,
403  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,
404  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,
405  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,
406  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,
407  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,
408  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,
409  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};
410 
411 
413 const unsigned char sAmbiguousFace[256] = {
414  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,
415  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,
416  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,
417  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,
418  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,
419  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,
420  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,
421  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};
422 
423 
427 const unsigned char sEdgeGroupTable[256][13] = {
428  {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},
429  {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},
430  {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},
431  {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},
432  {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},
433  {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},
434  {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},
435  {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},
436  {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},
437  {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},
438  {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},
439  {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},
440  {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},
441  {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},
442  {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},
443  {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},
444  {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},
445  {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},
446  {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},
447  {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},
448  {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},
449  {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},
450  {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},
451  {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},
452  {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},
453  {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},
454  {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},
455  {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},
456  {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},
457  {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},
458  {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},
459  {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},
460  {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},
461  {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},
462  {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},
463  {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},
464  {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},
465  {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},
466  {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},
467  {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},
468  {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},
469  {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},
470  {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},
471  {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},
472  {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},
473  {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},
474  {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},
475  {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},
476  {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},
477  {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},
478  {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},
479  {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},
480  {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},
481  {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},
482  {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},
483  {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},
484  {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},
485  {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},
486  {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},
487  {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},
488  {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},
489  {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},
490  {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},
491  {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},
492  {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},
493  {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},
494  {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},
495  {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},
496  {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},
497  {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},
498  {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},
499  {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},
500  {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},
501  {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},
502  {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},
503  {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},
504  {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},
505  {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},
506  {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},
507  {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},
508  {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},
509  {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},
510  {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},
511  {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},
512  {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},
513  {0,0,0,0,0,0,0,0,0,0,0,0,0}};
514 
515 
517 
518 inline bool
520  const Vec3d& p0, const Vec3d& p1,
521  const Vec3d& p2, const Vec3d& p3,
522  double epsilon = 0.001)
523 {
524  // compute representative plane
525  Vec3d normal = (p2-p0).cross(p1-p3);
526  normal.normalize();
527  const Vec3d centroid = (p0 + p1 + p2 + p3);
528  const double d = centroid.dot(normal) * 0.25;
529 
530 
531  // test vertice distance to plane
532  double absDist = std::abs(p0.dot(normal) - d);
533  if (absDist > epsilon) return false;
534 
535  absDist = std::abs(p1.dot(normal) - d);
536  if (absDist > epsilon) return false;
537 
538  absDist = std::abs(p2.dot(normal) - d);
539  if (absDist > epsilon) return false;
540 
541  absDist = std::abs(p3.dot(normal) - d);
542  if (absDist > epsilon) return false;
543 
544  return true;
545 }
546 
547 
549 
550 
553 
554 enum { MASK_FIRST_10_BITS = 0x000003FF, MASK_DIRTY_BIT = 0x80000000, MASK_INVALID_BIT = 0x40000000 };
555 
556 inline uint32_t
557 packPoint(const Vec3d& v)
558 {
559  uint32_t data = 0;
560 
561  // values are expected to be in the [0.0 to 1.0] range.
562  assert(!(v.x() > 1.0) && !(v.y() > 1.0) && !(v.z() > 1.0));
563  assert(!(v.x() < 0.0) && !(v.y() < 0.0) && !(v.z() < 0.0));
564 
565  data |= (uint32_t(v.x() * 1023.0) & MASK_FIRST_10_BITS) << 20;
566  data |= (uint32_t(v.y() * 1023.0) & MASK_FIRST_10_BITS) << 10;
567  data |= (uint32_t(v.z() * 1023.0) & MASK_FIRST_10_BITS);
568 
569  return data;
570 }
571 
572 inline Vec3d
573 unpackPoint(uint32_t data)
574 {
575  Vec3d v;
576  v.z() = double(data & MASK_FIRST_10_BITS) * 0.0009775171;
577  data = data >> 10;
578  v.y() = double(data & MASK_FIRST_10_BITS) * 0.0009775171;
579  data = data >> 10;
580  v.x() = double(data & MASK_FIRST_10_BITS) * 0.0009775171;
581 
582  return v;
583 }
584 
586 
588 
589 
592 template<typename AccessorT>
593 inline unsigned char
594 evalCellSigns(const AccessorT& accessor, const Coord& ijk, typename AccessorT::ValueType iso)
595 {
596  unsigned signs = 0;
597  Coord coord = ijk; // i, j, k
598  if (accessor.getValue(coord) < iso) signs |= 1u;
599  coord[0] += 1; // i+1, j, k
600  if (accessor.getValue(coord) < iso) signs |= 2u;
601  coord[2] += 1; // i+1, j, k+1
602  if (accessor.getValue(coord) < iso) signs |= 4u;
603  coord[0] = ijk[0]; // i, j, k+1
604  if (accessor.getValue(coord) < iso) signs |= 8u;
605  coord[1] += 1; coord[2] = ijk[2]; // i, j+1, k
606  if (accessor.getValue(coord) < iso) signs |= 16u;
607  coord[0] += 1; // i+1, j+1, k
608  if (accessor.getValue(coord) < iso) signs |= 32u;
609  coord[2] += 1; // i+1, j+1, k+1
610  if (accessor.getValue(coord) < iso) signs |= 64u;
611  coord[0] = ijk[0]; // i, j+1, k+1
612  if (accessor.getValue(coord) < iso) signs |= 128u;
613  return signs;
614 }
615 
616 
619 template<typename LeafT>
620 inline unsigned char
621 evalCellSigns(const LeafT& leaf, const Index offset, typename LeafT::ValueType iso)
622 {
623  unsigned char signs = 0;
624 
625  // i, j, k
626  if (leaf.getValue(offset) < iso) signs |= 1u;
627 
628  // i, j, k+1
629  if (leaf.getValue(offset + 1) < iso) signs |= 8u;
630 
631  // i, j+1, k
632  if (leaf.getValue(offset + LeafT::DIM) < iso) signs |= 16u;
633 
634  // i, j+1, k+1
635  if (leaf.getValue(offset + LeafT::DIM + 1) < iso) signs |= 128u;
636 
637  // i+1, j, k
638  if (leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) ) < iso) signs |= 2u;
639 
640  // i+1, j, k+1
641  if (leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + 1) < iso) signs |= 4u;
642 
643  // i+1, j+1, k
644  if (leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM) < iso) signs |= 32u;
645 
646  // i+1, j+1, k+1
647  if (leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM + 1) < iso) signs |= 64u;
648 
649  return signs;
650 }
651 
652 
655 template<class AccessorT>
656 inline void
657 correctCellSigns(unsigned char& signs, unsigned char face,
658  const AccessorT& acc, Coord ijk, typename AccessorT::ValueType iso)
659 {
660  if (face == 1) {
661  ijk[2] -= 1;
662  if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 3) signs = ~signs;
663  } else if (face == 3) {
664  ijk[2] += 1;
665  if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 1) signs = ~signs;
666  } else if (face == 2) {
667  ijk[0] += 1;
668  if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 4) signs = ~signs;
669  } else if (face == 4) {
670  ijk[0] -= 1;
671  if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 2) signs = ~signs;
672  } else if (face == 5) {
673  ijk[1] -= 1;
674  if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 6) signs = ~signs;
675  } else if (face == 6) {
676  ijk[1] += 1;
677  if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 5) signs = ~signs;
678  }
679 }
680 
681 
682 template<class AccessorT>
683 inline bool
684 isNonManifold(const AccessorT& accessor, const Coord& ijk,
685  typename AccessorT::ValueType isovalue, const int dim)
686 {
687  int hDim = dim >> 1;
688  bool m, p[8]; // Corner signs
689 
690  Coord coord = ijk; // i, j, k
691  p[0] = accessor.getValue(coord) < isovalue;
692  coord[0] += dim; // i+dim, j, k
693  p[1] = accessor.getValue(coord) < isovalue;
694  coord[2] += dim; // i+dim, j, k+dim
695  p[2] = accessor.getValue(coord) < isovalue;
696  coord[0] = ijk[0]; // i, j, k+dim
697  p[3] = accessor.getValue(coord) < isovalue;
698  coord[1] += dim; coord[2] = ijk[2]; // i, j+dim, k
699  p[4] = accessor.getValue(coord) < isovalue;
700  coord[0] += dim; // i+dim, j+dim, k
701  p[5] = accessor.getValue(coord) < isovalue;
702  coord[2] += dim; // i+dim, j+dim, k+dim
703  p[6] = accessor.getValue(coord) < isovalue;
704  coord[0] = ijk[0]; // i, j+dim, k+dim
705  p[7] = accessor.getValue(coord) < isovalue;
706 
707  // Check if the corner sign configuration is ambiguous
708  unsigned signs = 0;
709  if (p[0]) signs |= 1u;
710  if (p[1]) signs |= 2u;
711  if (p[2]) signs |= 4u;
712  if (p[3]) signs |= 8u;
713  if (p[4]) signs |= 16u;
714  if (p[5]) signs |= 32u;
715  if (p[6]) signs |= 64u;
716  if (p[7]) signs |= 128u;
717  if (!sAdaptable[signs]) return true;
718 
719  // Manifold check
720 
721  // Evaluate edges
722  int i = ijk[0], ip = ijk[0] + hDim, ipp = ijk[0] + dim;
723  int j = ijk[1], jp = ijk[1] + hDim, jpp = ijk[1] + dim;
724  int k = ijk[2], kp = ijk[2] + hDim, kpp = ijk[2] + dim;
725 
726  // edge 1
727  coord.reset(ip, j, k);
728  m = accessor.getValue(coord) < isovalue;
729  if (p[0] != m && p[1] != m) return true;
730 
731  // edge 2
732  coord.reset(ipp, j, kp);
733  m = accessor.getValue(coord) < isovalue;
734  if (p[1] != m && p[2] != m) return true;
735 
736  // edge 3
737  coord.reset(ip, j, kpp);
738  m = accessor.getValue(coord) < isovalue;
739  if (p[2] != m && p[3] != m) return true;
740 
741  // edge 4
742  coord.reset(i, j, kp);
743  m = accessor.getValue(coord) < isovalue;
744  if (p[0] != m && p[3] != m) return true;
745 
746  // edge 5
747  coord.reset(ip, jpp, k);
748  m = accessor.getValue(coord) < isovalue;
749  if (p[4] != m && p[5] != m) return true;
750 
751  // edge 6
752  coord.reset(ipp, jpp, kp);
753  m = accessor.getValue(coord) < isovalue;
754  if (p[5] != m && p[6] != m) return true;
755 
756  // edge 7
757  coord.reset(ip, jpp, kpp);
758  m = accessor.getValue(coord) < isovalue;
759  if (p[6] != m && p[7] != m) return true;
760 
761  // edge 8
762  coord.reset(i, jpp, kp);
763  m = accessor.getValue(coord) < isovalue;
764  if (p[7] != m && p[4] != m) return true;
765 
766  // edge 9
767  coord.reset(i, jp, k);
768  m = accessor.getValue(coord) < isovalue;
769  if (p[0] != m && p[4] != m) return true;
770 
771  // edge 10
772  coord.reset(ipp, jp, k);
773  m = accessor.getValue(coord) < isovalue;
774  if (p[1] != m && p[5] != m) return true;
775 
776  // edge 11
777  coord.reset(ipp, jp, kpp);
778  m = accessor.getValue(coord) < isovalue;
779  if (p[2] != m && p[6] != m) return true;
780 
781 
782  // edge 12
783  coord.reset(i, jp, kpp);
784  m = accessor.getValue(coord) < isovalue;
785  if (p[3] != m && p[7] != m) return true;
786 
787 
788  // Evaluate faces
789 
790  // face 1
791  coord.reset(ip, jp, k);
792  m = accessor.getValue(coord) < isovalue;
793  if (p[0] != m && p[1] != m && p[4] != m && p[5] != m) return true;
794 
795  // face 2
796  coord.reset(ipp, jp, kp);
797  m = accessor.getValue(coord) < isovalue;
798  if (p[1] != m && p[2] != m && p[5] != m && p[6] != m) return true;
799 
800  // face 3
801  coord.reset(ip, jp, kpp);
802  m = accessor.getValue(coord) < isovalue;
803  if (p[2] != m && p[3] != m && p[6] != m && p[7] != m) return true;
804 
805  // face 4
806  coord.reset(i, jp, kp);
807  m = accessor.getValue(coord) < isovalue;
808  if (p[0] != m && p[3] != m && p[4] != m && p[7] != m) return true;
809 
810  // face 5
811  coord.reset(ip, j, kp);
812  m = accessor.getValue(coord) < isovalue;
813  if (p[0] != m && p[1] != m && p[2] != m && p[3] != m) return true;
814 
815  // face 6
816  coord.reset(ip, jpp, kp);
817  m = accessor.getValue(coord) < isovalue;
818  if (p[4] != m && p[5] != m && p[6] != m && p[7] != m) return true;
819 
820  // test cube center
821  coord.reset(ip, jp, kp);
822  m = accessor.getValue(coord) < isovalue;
823  if (p[0] != m && p[1] != m && p[2] != m && p[3] != m &&
824  p[4] != m && p[5] != m && p[6] != m && p[7] != m) return true;
825 
826  return false;
827 }
828 
829 
831 
832 
833 template <class LeafType>
834 inline void
835 mergeVoxels(LeafType& leaf, const Coord& start, int dim, int regionId)
836 {
837  Coord ijk, end = start;
838  end[0] += dim;
839  end[1] += dim;
840  end[2] += dim;
841 
842  for (ijk[0] = start[0]; ijk[0] < end[0]; ++ijk[0]) {
843  for (ijk[1] = start[1]; ijk[1] < end[1]; ++ijk[1]) {
844  for (ijk[2] = start[2]; ijk[2] < end[2]; ++ijk[2]) {
845  if(leaf.isValueOn(ijk)) leaf.setValue(ijk, regionId);
846  }
847  }
848  }
849 }
850 
851 
852 // Note that we must use ValueType::value_type or else Visual C++ gets confused
853 // thinking that it is a constructor.
854 template <class LeafType>
855 inline bool
856 isMergable(LeafType& leaf, const Coord& start, int dim,
857  typename LeafType::ValueType::value_type adaptivity)
858 {
859  if (adaptivity < 1e-6) return false;
860 
861  typedef typename LeafType::ValueType VecT;
862  Coord ijk, end = start;
863  end[0] += dim;
864  end[1] += dim;
865  end[2] += dim;
866 
867  std::vector<VecT> norms;
868  for (ijk[0] = start[0]; ijk[0] < end[0]; ++ijk[0]) {
869  for (ijk[1] = start[1]; ijk[1] < end[1]; ++ijk[1]) {
870  for (ijk[2] = start[2]; ijk[2] < end[2]; ++ijk[2]) {
871 
872  if(!leaf.isValueOn(ijk)) continue;
873  norms.push_back(leaf.getValue(ijk));
874  }
875  }
876  }
877 
878  size_t N = norms.size();
879  for (size_t ni = 0; ni < N; ++ni) {
880  VecT n_i = norms[ni];
881  for (size_t nj = 0; nj < N; ++nj) {
882  VecT n_j = norms[nj];
883  if ((1.0 - n_i.dot(n_j)) > adaptivity) return false;
884  }
885  }
886  return true;
887 }
888 
889 
891 
892 
893 template<typename TreeT, typename LeafManagerT>
894 class SignData
895 {
896 public:
897  typedef typename TreeT::ValueType ValueT;
899 
900  typedef typename TreeT::template ValueConverter<int>::Type IntTreeT;
902 
903  typedef typename TreeT::template ValueConverter<Int16>::Type Int16TreeT;
905 
907 
908 
909  SignData(const TreeT& distTree, const LeafManagerT& leafs, ValueT iso);
910 
911  void run(bool threaded = true);
912 
913  typename Int16TreeT::Ptr signTree() const { return mSignTree; }
914  typename IntTreeT::Ptr idxTree() const { return mIdxTree; }
915 
917 
918  SignData(SignData&, tbb::split);
919  void operator()(const tbb::blocked_range<size_t>&);
920  void join(const SignData& rhs)
921  {
922  mSignTree->merge(*rhs.mSignTree);
923  mIdxTree->merge(*rhs.mIdxTree);
924  }
925 
926 private:
927 
928  const TreeT& mDistTree;
929  AccessorT mDistAcc;
930 
931  const LeafManagerT& mLeafs;
932  ValueT mIsovalue;
933 
934  typename Int16TreeT::Ptr mSignTree;
935  Int16AccessorT mSignAcc;
936 
937  typename IntTreeT::Ptr mIdxTree;
938  IntAccessorT mIdxAcc;
939 
940 };
941 
942 
943 template<typename TreeT, typename LeafManagerT>
945  const LeafManagerT& leafs, ValueT iso)
946  : mDistTree(distTree)
947  , mDistAcc(mDistTree)
948  , mLeafs(leafs)
949  , mIsovalue(iso)
950  , mSignTree(new Int16TreeT(0))
951  , mSignAcc(*mSignTree)
952  , mIdxTree(new IntTreeT(int(util::INVALID_IDX)))
953  , mIdxAcc(*mIdxTree)
954 {
955 }
956 
957 
958 template<typename TreeT, typename LeafManagerT>
960  : mDistTree(rhs.mDistTree)
961  , mDistAcc(mDistTree)
962  , mLeafs(rhs.mLeafs)
963  , mIsovalue(rhs.mIsovalue)
964  , mSignTree(new Int16TreeT(0))
965  , mSignAcc(*mSignTree)
966  , mIdxTree(new IntTreeT(int(util::INVALID_IDX)))
967  , mIdxAcc(*mIdxTree)
968 {
969 }
970 
971 
972 template<typename TreeT, typename LeafManagerT>
973 void
975 {
976  if (threaded) tbb::parallel_reduce(mLeafs.getRange(), *this);
977  else (*this)(mLeafs.getRange());
978 }
979 
980 template<typename TreeT, typename LeafManagerT>
981 void
982 SignData<TreeT, LeafManagerT>::operator()(const tbb::blocked_range<size_t>& range)
983 {
984  typedef typename Int16TreeT::LeafNodeType Int16LeafT;
985  typedef typename IntTreeT::LeafNodeType IntLeafT;
986  typename LeafManagerT::TreeType::LeafNodeType::ValueOnCIter iter;
987  unsigned char signs, face;
988  Coord ijk, coord;
989 
990  typename internal::UniquePtr<Int16LeafT>::type signLeafPt(new Int16LeafT(ijk, 0));
991 
992  for (size_t n = range.begin(); n != range.end(); ++n) {
993 
994  bool collectedData = false;
995 
996  coord = mLeafs.leaf(n).origin();
997 
998  if (!signLeafPt.get()) signLeafPt.reset(new Int16LeafT(coord, 0));
999  else signLeafPt->setOrigin(coord);
1000 
1001  const typename TreeT::LeafNodeType *leafPt = mDistAcc.probeConstLeaf(coord);
1002 
1003  coord.offset(TreeT::LeafNodeType::DIM - 1);
1004 
1005  for (iter = mLeafs.leaf(n).cbeginValueOn(); iter; ++iter) {
1006 
1007  ijk = iter.getCoord();
1008 
1009  if (leafPt && ijk[0] < coord[0] && ijk[1] < coord[1] && ijk[2] < coord[2]) {
1010  signs = evalCellSigns(*leafPt, iter.pos(), mIsovalue);
1011  } else {
1012  signs = evalCellSigns(mDistAcc, ijk, mIsovalue);
1013  }
1014 
1015  if (signs != 0 && signs != 0xFF) {
1016  Int16 flags = (signs & 0x1) ? INSIDE : 0;
1017 
1018  if (bool(signs & 0x1) != bool(signs & 0x2)) flags |= XEDGE;
1019  if (bool(signs & 0x1) != bool(signs & 0x10)) flags |= YEDGE;
1020  if (bool(signs & 0x1) != bool(signs & 0x8)) flags |= ZEDGE;
1021 
1022  face = internal::sAmbiguousFace[signs];
1023  if (face != 0) correctCellSigns(signs, face, mDistAcc, ijk, mIsovalue);
1024 
1025  flags |= Int16(signs);
1026 
1027  signLeafPt->setValue(ijk, flags);
1028  collectedData = true;
1029  }
1030  }
1031 
1032  if (collectedData) {
1033 
1034  IntLeafT* idxLeaf = mIdxAcc.touchLeaf(coord);
1035  idxLeaf->topologyUnion(*signLeafPt);
1036  typename IntLeafT::ValueOnIter it = idxLeaf->beginValueOn();
1037  for (; it; ++it) {
1038  it.setValue(0);
1039  }
1040 
1041  mSignAcc.addLeaf(signLeafPt.release());
1042  }
1043  }
1044 }
1045 
1046 
1048 
1049 
1052 {
1053 public:
1054  CountPoints(std::vector<size_t>& pointList) : mPointList(pointList) {}
1055 
1056  template <typename LeafNodeType>
1057  void operator()(LeafNodeType &leaf, size_t leafIndex) const
1058  {
1059  size_t points = 0;
1060 
1061  typename LeafNodeType::ValueOnCIter iter = leaf.cbeginValueOn();
1062  for (; iter; ++iter) {
1063  points += size_t(sEdgeGroupTable[(SIGNS & iter.getValue())][0]);
1064  }
1065 
1066  mPointList[leafIndex] = points;
1067  }
1068 
1069 private:
1070  std::vector<size_t>& mPointList;
1071 };
1072 
1073 
1075 template<typename Int16TreeT>
1077 {
1078 public:
1080 
1081  MapPoints(std::vector<size_t>& pointList, const Int16TreeT& signTree)
1082  : mPointList(pointList)
1083  , mSignAcc(signTree)
1084  {
1085  }
1086 
1087  template <typename LeafNodeType>
1088  void operator()(LeafNodeType &leaf, size_t leafIndex) const
1089  {
1090  size_t ptnIdx = mPointList[leafIndex];
1091  typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
1092 
1093  const typename Int16TreeT::LeafNodeType *signLeafPt =
1094  mSignAcc.probeConstLeaf(leaf.origin());
1095 
1096  for (; iter; ++iter) {
1097  iter.setValue(ptnIdx);
1098  unsigned signs = SIGNS & signLeafPt->getValue(iter.pos());
1099  ptnIdx += size_t(sEdgeGroupTable[signs][0]);
1100  }
1101  }
1102 
1103 private:
1104  std::vector<size_t>& mPointList;
1105  Int16AccessorT mSignAcc;
1106 };
1107 
1108 
1110 template<typename IntTreeT>
1112 {
1113 public:
1115  typedef typename IntTreeT::LeafNodeType IntLeafT;
1116 
1117  CountRegions(IntTreeT& idxTree, std::vector<size_t>& regions)
1118  : mIdxAcc(idxTree)
1119  , mRegions(regions)
1120  {
1121  }
1122 
1123  template <typename LeafNodeType>
1124  void operator()(LeafNodeType &leaf, size_t leafIndex) const
1125  {
1126 
1127  size_t regions = 0;
1128 
1129  IntLeafT tmpLeaf(*mIdxAcc.probeConstLeaf(leaf.origin()));
1130 
1131  typename IntLeafT::ValueOnIter iter = tmpLeaf.beginValueOn();
1132  for (; iter; ++iter) {
1133  if(iter.getValue() == 0) {
1134  iter.setValueOff();
1135  regions += size_t(sEdgeGroupTable[(SIGNS & leaf.getValue(iter.pos()))][0]);
1136  }
1137  }
1138 
1139  int onVoxelCount = int(tmpLeaf.onVoxelCount());
1140  while (onVoxelCount > 0) {
1141  ++regions;
1142  iter = tmpLeaf.beginValueOn();
1143  int regionId = iter.getValue();
1144  for (; iter; ++iter) {
1145  if (iter.getValue() == regionId) {
1146  iter.setValueOff();
1147  --onVoxelCount;
1148  }
1149  }
1150  }
1151 
1152  mRegions[leafIndex] = regions;
1153  }
1154 
1155 private:
1156  IntAccessorT mIdxAcc;
1157  std::vector<size_t>& mRegions;
1158 };
1159 
1160 
1162 
1163 
1164 // @brief linear interpolation.
1165 inline double evalRoot(double v0, double v1, double iso) { return (iso - v0) / (v1 - v0); }
1166 
1167 
1169 template<typename LeafT>
1170 inline void
1171 collectCornerValues(const LeafT& leaf, const Index offset, std::vector<double>& values)
1172 {
1173  values[0] = double(leaf.getValue(offset)); // i, j, k
1174  values[3] = double(leaf.getValue(offset + 1)); // i, j, k+1
1175  values[4] = double(leaf.getValue(offset + LeafT::DIM)); // i, j+1, k
1176  values[7] = double(leaf.getValue(offset + LeafT::DIM + 1)); // i, j+1, k+1
1177  values[1] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM))); // i+1, j, k
1178  values[2] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + 1)); // i+1, j, k+1
1179  values[5] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM)); // i+1, j+1, k
1180  values[6] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM + 1)); // i+1, j+1, k+1
1181 }
1182 
1183 
1185 template<typename AccessorT>
1186 inline void
1187 collectCornerValues(const AccessorT& acc, const Coord& ijk, std::vector<double>& values)
1188 {
1189  Coord coord = ijk;
1190  values[0] = double(acc.getValue(coord)); // i, j, k
1191 
1192  coord[0] += 1;
1193  values[1] = double(acc.getValue(coord)); // i+1, j, k
1194 
1195  coord[2] += 1;
1196  values[2] = double(acc.getValue(coord)); // i+i, j, k+1
1197 
1198  coord[0] = ijk[0];
1199  values[3] = double(acc.getValue(coord)); // i, j, k+1
1200 
1201  coord[1] += 1; coord[2] = ijk[2];
1202  values[4] = double(acc.getValue(coord)); // i, j+1, k
1203 
1204  coord[0] += 1;
1205  values[5] = double(acc.getValue(coord)); // i+1, j+1, k
1206 
1207  coord[2] += 1;
1208  values[6] = double(acc.getValue(coord)); // i+1, j+1, k+1
1209 
1210  coord[0] = ijk[0];
1211  values[7] = double(acc.getValue(coord)); // i, j+1, k+1
1212 }
1213 
1214 
1216 inline Vec3d
1217 computePoint(const std::vector<double>& values, unsigned char signs,
1218  unsigned char edgeGroup, double iso)
1219 {
1220  Vec3d avg(0.0, 0.0, 0.0);
1221  int samples = 0;
1222 
1223  if (sEdgeGroupTable[signs][1] == edgeGroup) { // Edged: 0 - 1
1224  avg[0] += evalRoot(values[0], values[1], iso);
1225  ++samples;
1226  }
1227 
1228  if (sEdgeGroupTable[signs][2] == edgeGroup) { // Edged: 1 - 2
1229  avg[0] += 1.0;
1230  avg[2] += evalRoot(values[1], values[2], iso);
1231  ++samples;
1232  }
1233 
1234  if (sEdgeGroupTable[signs][3] == edgeGroup) { // Edged: 3 - 2
1235  avg[0] += evalRoot(values[3], values[2], iso);
1236  avg[2] += 1.0;
1237  ++samples;
1238  }
1239 
1240  if (sEdgeGroupTable[signs][4] == edgeGroup) { // Edged: 0 - 3
1241  avg[2] += evalRoot(values[0], values[3], iso);
1242  ++samples;
1243  }
1244 
1245  if (sEdgeGroupTable[signs][5] == edgeGroup) { // Edged: 4 - 5
1246  avg[0] += evalRoot(values[4], values[5], iso);
1247  avg[1] += 1.0;
1248  ++samples;
1249  }
1250 
1251  if (sEdgeGroupTable[signs][6] == edgeGroup) { // Edged: 5 - 6
1252  avg[0] += 1.0;
1253  avg[1] += 1.0;
1254  avg[2] += evalRoot(values[5], values[6], iso);
1255  ++samples;
1256  }
1257 
1258  if (sEdgeGroupTable[signs][7] == edgeGroup) { // Edged: 7 - 6
1259  avg[0] += evalRoot(values[7], values[6], iso);
1260  avg[1] += 1.0;
1261  avg[2] += 1.0;
1262  ++samples;
1263  }
1264 
1265  if (sEdgeGroupTable[signs][8] == edgeGroup) { // Edged: 4 - 7
1266  avg[1] += 1.0;
1267  avg[2] += evalRoot(values[4], values[7], iso);
1268  ++samples;
1269  }
1270 
1271  if (sEdgeGroupTable[signs][9] == edgeGroup) { // Edged: 0 - 4
1272  avg[1] += evalRoot(values[0], values[4], iso);
1273  ++samples;
1274  }
1275 
1276  if (sEdgeGroupTable[signs][10] == edgeGroup) { // Edged: 1 - 5
1277  avg[0] += 1.0;
1278  avg[1] += evalRoot(values[1], values[5], iso);
1279  ++samples;
1280  }
1281 
1282  if (sEdgeGroupTable[signs][11] == edgeGroup) { // Edged: 2 - 6
1283  avg[0] += 1.0;
1284  avg[1] += evalRoot(values[2], values[6], iso);
1285  avg[2] += 1.0;
1286  ++samples;
1287  }
1288 
1289  if (sEdgeGroupTable[signs][12] == edgeGroup) { // Edged: 3 - 7
1290  avg[1] += evalRoot(values[3], values[7], iso);
1291  avg[2] += 1.0;
1292  ++samples;
1293  }
1294 
1295  if (samples > 1) {
1296  double w = 1.0 / double(samples);
1297  avg[0] *= w;
1298  avg[1] *= w;
1299  avg[2] *= w;
1300  }
1301 
1302  return avg;
1303 }
1304 
1305 
1308 inline int
1309 computeMaskedPoint(Vec3d& avg, const std::vector<double>& values, unsigned char signs,
1310  unsigned char signsMask, unsigned char edgeGroup, double iso)
1311 {
1312  avg = Vec3d(0.0, 0.0, 0.0);
1313  int samples = 0;
1314 
1315  if (sEdgeGroupTable[signs][1] == edgeGroup
1316  && sEdgeGroupTable[signsMask][1] == 0) { // Edged: 0 - 1
1317  avg[0] += evalRoot(values[0], values[1], iso);
1318  ++samples;
1319  }
1320 
1321  if (sEdgeGroupTable[signs][2] == edgeGroup
1322  && sEdgeGroupTable[signsMask][2] == 0) { // Edged: 1 - 2
1323  avg[0] += 1.0;
1324  avg[2] += evalRoot(values[1], values[2], iso);
1325  ++samples;
1326  }
1327 
1328  if (sEdgeGroupTable[signs][3] == edgeGroup
1329  && sEdgeGroupTable[signsMask][3] == 0) { // Edged: 3 - 2
1330  avg[0] += evalRoot(values[3], values[2], iso);
1331  avg[2] += 1.0;
1332  ++samples;
1333  }
1334 
1335  if (sEdgeGroupTable[signs][4] == edgeGroup
1336  && sEdgeGroupTable[signsMask][4] == 0) { // Edged: 0 - 3
1337  avg[2] += evalRoot(values[0], values[3], iso);
1338  ++samples;
1339  }
1340 
1341  if (sEdgeGroupTable[signs][5] == edgeGroup
1342  && sEdgeGroupTable[signsMask][5] == 0) { // Edged: 4 - 5
1343  avg[0] += evalRoot(values[4], values[5], iso);
1344  avg[1] += 1.0;
1345  ++samples;
1346  }
1347 
1348  if (sEdgeGroupTable[signs][6] == edgeGroup
1349  && sEdgeGroupTable[signsMask][6] == 0) { // Edged: 5 - 6
1350  avg[0] += 1.0;
1351  avg[1] += 1.0;
1352  avg[2] += evalRoot(values[5], values[6], iso);
1353  ++samples;
1354  }
1355 
1356  if (sEdgeGroupTable[signs][7] == edgeGroup
1357  && sEdgeGroupTable[signsMask][7] == 0) { // Edged: 7 - 6
1358  avg[0] += evalRoot(values[7], values[6], iso);
1359  avg[1] += 1.0;
1360  avg[2] += 1.0;
1361  ++samples;
1362  }
1363 
1364  if (sEdgeGroupTable[signs][8] == edgeGroup
1365  && sEdgeGroupTable[signsMask][8] == 0) { // Edged: 4 - 7
1366  avg[1] += 1.0;
1367  avg[2] += evalRoot(values[4], values[7], iso);
1368  ++samples;
1369  }
1370 
1371  if (sEdgeGroupTable[signs][9] == edgeGroup
1372  && sEdgeGroupTable[signsMask][9] == 0) { // Edged: 0 - 4
1373  avg[1] += evalRoot(values[0], values[4], iso);
1374  ++samples;
1375  }
1376 
1377  if (sEdgeGroupTable[signs][10] == edgeGroup
1378  && sEdgeGroupTable[signsMask][10] == 0) { // Edged: 1 - 5
1379  avg[0] += 1.0;
1380  avg[1] += evalRoot(values[1], values[5], iso);
1381  ++samples;
1382  }
1383 
1384  if (sEdgeGroupTable[signs][11] == edgeGroup
1385  && sEdgeGroupTable[signsMask][11] == 0) { // Edged: 2 - 6
1386  avg[0] += 1.0;
1387  avg[1] += evalRoot(values[2], values[6], iso);
1388  avg[2] += 1.0;
1389  ++samples;
1390  }
1391 
1392  if (sEdgeGroupTable[signs][12] == edgeGroup
1393  && sEdgeGroupTable[signsMask][12] == 0) { // Edged: 3 - 7
1394  avg[1] += evalRoot(values[3], values[7], iso);
1395  avg[2] += 1.0;
1396  ++samples;
1397  }
1398 
1399  if (samples > 1) {
1400  double w = 1.0 / double(samples);
1401  avg[0] *= w;
1402  avg[1] *= w;
1403  avg[2] *= w;
1404  }
1405 
1406  return samples;
1407 }
1408 
1409 
1412 inline Vec3d
1413 computeWeightedPoint(const Vec3d& p, const std::vector<double>& values,
1414  unsigned char signs, unsigned char edgeGroup, double iso)
1415 {
1416  std::vector<Vec3d> samples;
1417  samples.reserve(8);
1418 
1419  std::vector<double> weights;
1420  weights.reserve(8);
1421 
1422  Vec3d avg(0.0, 0.0, 0.0);
1423 
1424  if (sEdgeGroupTable[signs][1] == edgeGroup) { // Edged: 0 - 1
1425  avg[0] = evalRoot(values[0], values[1], iso);
1426  avg[1] = 0.0;
1427  avg[2] = 0.0;
1428 
1429  samples.push_back(avg);
1430  weights.push_back((avg-p).lengthSqr());
1431  }
1432 
1433  if (sEdgeGroupTable[signs][2] == edgeGroup) { // Edged: 1 - 2
1434  avg[0] = 1.0;
1435  avg[1] = 0.0;
1436  avg[2] = evalRoot(values[1], values[2], iso);
1437 
1438  samples.push_back(avg);
1439  weights.push_back((avg-p).lengthSqr());
1440  }
1441 
1442  if (sEdgeGroupTable[signs][3] == edgeGroup) { // Edged: 3 - 2
1443  avg[0] = evalRoot(values[3], values[2], iso);
1444  avg[1] = 0.0;
1445  avg[2] = 1.0;
1446 
1447  samples.push_back(avg);
1448  weights.push_back((avg-p).lengthSqr());
1449  }
1450 
1451  if (sEdgeGroupTable[signs][4] == edgeGroup) { // Edged: 0 - 3
1452  avg[0] = 0.0;
1453  avg[1] = 0.0;
1454  avg[2] = evalRoot(values[0], values[3], iso);
1455 
1456  samples.push_back(avg);
1457  weights.push_back((avg-p).lengthSqr());
1458  }
1459 
1460  if (sEdgeGroupTable[signs][5] == edgeGroup) { // Edged: 4 - 5
1461  avg[0] = evalRoot(values[4], values[5], iso);
1462  avg[1] = 1.0;
1463  avg[2] = 0.0;
1464 
1465  samples.push_back(avg);
1466  weights.push_back((avg-p).lengthSqr());
1467  }
1468 
1469  if (sEdgeGroupTable[signs][6] == edgeGroup) { // Edged: 5 - 6
1470  avg[0] = 1.0;
1471  avg[1] = 1.0;
1472  avg[2] = evalRoot(values[5], values[6], iso);
1473 
1474  samples.push_back(avg);
1475  weights.push_back((avg-p).lengthSqr());
1476  }
1477 
1478  if (sEdgeGroupTable[signs][7] == edgeGroup) { // Edged: 7 - 6
1479  avg[0] = evalRoot(values[7], values[6], iso);
1480  avg[1] = 1.0;
1481  avg[2] = 1.0;
1482 
1483  samples.push_back(avg);
1484  weights.push_back((avg-p).lengthSqr());
1485  }
1486 
1487  if (sEdgeGroupTable[signs][8] == edgeGroup) { // Edged: 4 - 7
1488  avg[0] = 0.0;
1489  avg[1] = 1.0;
1490  avg[2] = evalRoot(values[4], values[7], iso);
1491 
1492  samples.push_back(avg);
1493  weights.push_back((avg-p).lengthSqr());
1494  }
1495 
1496  if (sEdgeGroupTable[signs][9] == edgeGroup) { // Edged: 0 - 4
1497  avg[0] = 0.0;
1498  avg[1] = evalRoot(values[0], values[4], iso);
1499  avg[2] = 0.0;
1500 
1501  samples.push_back(avg);
1502  weights.push_back((avg-p).lengthSqr());
1503  }
1504 
1505  if (sEdgeGroupTable[signs][10] == edgeGroup) { // Edged: 1 - 5
1506  avg[0] = 1.0;
1507  avg[1] = evalRoot(values[1], values[5], iso);
1508  avg[2] = 0.0;
1509 
1510  samples.push_back(avg);
1511  weights.push_back((avg-p).lengthSqr());
1512  }
1513 
1514  if (sEdgeGroupTable[signs][11] == edgeGroup) { // Edged: 2 - 6
1515  avg[0] = 1.0;
1516  avg[1] = evalRoot(values[2], values[6], iso);
1517  avg[2] = 1.0;
1518 
1519  samples.push_back(avg);
1520  weights.push_back((avg-p).lengthSqr());
1521  }
1522 
1523  if (sEdgeGroupTable[signs][12] == edgeGroup) { // Edged: 3 - 7
1524  avg[0] = 0.0;
1525  avg[1] = evalRoot(values[3], values[7], iso);
1526  avg[2] = 1.0;
1527 
1528  samples.push_back(avg);
1529  weights.push_back((avg-p).lengthSqr());
1530  }
1531 
1532 
1533  double minWeight = std::numeric_limits<double>::max();
1534  double maxWeight = -std::numeric_limits<double>::max();
1535 
1536  for (size_t i = 0, I = weights.size(); i < I; ++i) {
1537  minWeight = std::min(minWeight, weights[i]);
1538  maxWeight = std::max(maxWeight, weights[i]);
1539  }
1540 
1541  const double offset = maxWeight + minWeight * 0.1;
1542  for (size_t i = 0, I = weights.size(); i < I; ++i) {
1543  weights[i] = offset - weights[i];
1544  }
1545 
1546 
1547  double weightSum = 0.0;
1548  for (size_t i = 0, I = weights.size(); i < I; ++i) {
1549  weightSum += weights[i];
1550  }
1551 
1552  avg[0] = 0.0;
1553  avg[1] = 0.0;
1554  avg[2] = 0.0;
1555 
1556  if (samples.size() > 1) {
1557  for (size_t i = 0, I = samples.size(); i < I; ++i) {
1558  avg += samples[i] * (weights[i] / weightSum);
1559  }
1560  } else {
1561  avg = samples.front();
1562  }
1563 
1564  return avg;
1565 }
1566 
1567 
1570 inline void
1571 computeCellPoints(std::vector<Vec3d>& points,
1572  const std::vector<double>& values, unsigned char signs, double iso)
1573 {
1574  for (size_t n = 1, N = sEdgeGroupTable[signs][0] + 1; n < N; ++n) {
1575  points.push_back(computePoint(values, signs, n, iso));
1576  }
1577 }
1578 
1579 
1583 inline int
1584 matchEdgeGroup(unsigned char groupId, unsigned char lhsSigns, unsigned char rhsSigns)
1585 {
1586  int id = -1;
1587  for (size_t i = 1; i <= 12; ++i) {
1588  if (sEdgeGroupTable[lhsSigns][i] == groupId && sEdgeGroupTable[rhsSigns][i] != 0) {
1589  id = sEdgeGroupTable[rhsSigns][i];
1590  break;
1591  }
1592  }
1593  return id;
1594 }
1595 
1596 
1601 inline void
1602 computeCellPoints(std::vector<Vec3d>& points, std::vector<bool>& weightedPointMask,
1603  const std::vector<double>& lhsValues, const std::vector<double>& rhsValues,
1604  unsigned char lhsSigns, unsigned char rhsSigns,
1605  double iso, size_t pointIdx, const boost::scoped_array<uint32_t>& seamPoints)
1606 {
1607  for (size_t n = 1, N = sEdgeGroupTable[lhsSigns][0] + 1; n < N; ++n) {
1608 
1609  int id = matchEdgeGroup(n, lhsSigns, rhsSigns);
1610 
1611  if (id != -1) {
1612 
1613  const unsigned char e(id);
1614  uint32_t& quantizedPoint = seamPoints[pointIdx + (id - 1)];
1615 
1616  if ((quantizedPoint & MASK_DIRTY_BIT) && !(quantizedPoint & MASK_INVALID_BIT)) {
1617  Vec3d p = unpackPoint(quantizedPoint);
1618  points.push_back(computeWeightedPoint(p, rhsValues, rhsSigns, e, iso));
1619  weightedPointMask.push_back(true);
1620  } else {
1621  points.push_back(computePoint(rhsValues, rhsSigns, e, iso));
1622  weightedPointMask.push_back(false);
1623  }
1624 
1625  } else {
1626  points.push_back(computePoint(lhsValues, lhsSigns, n, iso));
1627  weightedPointMask.push_back(false);
1628  }
1629  }
1630 }
1631 
1632 
1633 template <typename TreeT, typename LeafManagerT>
1635 {
1636 public:
1638 
1639  typedef typename TreeT::template ValueConverter<int>::Type IntTreeT;
1642 
1643  typedef typename TreeT::template ValueConverter<Int16>::Type Int16TreeT;
1645 
1646  typedef boost::scoped_array<uint32_t> QuantizedPointList;
1647 
1649 
1650 
1651  GenPoints(const LeafManagerT& signLeafs, const TreeT& distTree,
1652  IntTreeT& idxTree, PointList& points, std::vector<size_t>& indices,
1653  const math::Transform& xform, double iso);
1654 
1655  void run(bool threaded = true);
1656 
1657  void setRefData(const Int16TreeT* refSignTree = NULL, const TreeT* refDistTree = NULL,
1658  IntTreeT* refIdxTree = NULL, const QuantizedPointList* seamPoints = NULL,
1659  std::vector<unsigned char>* mSeamPointMaskPt = NULL);
1660 
1662 
1663 
1664  void operator()(const tbb::blocked_range<size_t>&) const;
1665 
1666 private:
1667  const LeafManagerT& mSignLeafs;
1668 
1669  AccessorT mDistAcc;
1670  IntTreeT& mIdxTree;
1671 
1672  PointList& mPoints;
1673  std::vector<size_t>& mIndices;
1674  const math::Transform& mTransform;
1675  const double mIsovalue;
1676 
1677  // reference data
1678  const Int16TreeT *mRefSignTreePt;
1679  const TreeT* mRefDistTreePt;
1680  const IntTreeT* mRefIdxTreePt;
1681  const QuantizedPointList* mSeamPointsPt;
1682  std::vector<unsigned char>* mSeamPointMaskPt;
1683 };
1684 
1685 
1686 template <typename TreeT, typename LeafManagerT>
1687 GenPoints<TreeT, LeafManagerT>::GenPoints(const LeafManagerT& signLeafs,
1688  const TreeT& distTree, IntTreeT& idxTree, PointList& points,
1689  std::vector<size_t>& indices, const math::Transform& xform, double iso)
1690  : mSignLeafs(signLeafs)
1691  , mDistAcc(distTree)
1692  , mIdxTree(idxTree)
1693  , mPoints(points)
1694  , mIndices(indices)
1695  , mTransform(xform)
1696  , mIsovalue(iso)
1697  , mRefSignTreePt(NULL)
1698  , mRefDistTreePt(NULL)
1699  , mRefIdxTreePt(NULL)
1700  , mSeamPointsPt(NULL)
1701  , mSeamPointMaskPt(NULL)
1702 {
1703 }
1704 
1705 
1706 template <typename TreeT, typename LeafManagerT>
1707 void
1709 {
1710  if (threaded) tbb::parallel_for(mSignLeafs.getRange(), *this);
1711  else (*this)(mSignLeafs.getRange());
1712 }
1713 
1714 
1715 template <typename TreeT, typename LeafManagerT>
1716 void
1717 GenPoints<TreeT, LeafManagerT>::setRefData(const Int16TreeT *refSignTree, const TreeT *refDistTree,
1718  IntTreeT* refIdxTree, const QuantizedPointList* seamPoints, std::vector<unsigned char>* seamPointMask)
1719 {
1720  mRefSignTreePt = refSignTree;
1721  mRefDistTreePt = refDistTree;
1722  mRefIdxTreePt = refIdxTree;
1723  mSeamPointsPt = seamPoints;
1724  mSeamPointMaskPt = seamPointMask;
1725 }
1726 
1727 
1728 template <typename TreeT, typename LeafManagerT>
1729 void
1731  const tbb::blocked_range<size_t>& range) const
1732 {
1733  typename IntTreeT::LeafNodeType::ValueOnIter iter;
1734  unsigned char signs, refSigns;
1735  Index offset;
1736  Coord ijk, coord;
1737  std::vector<Vec3d> points(4);
1738  std::vector<bool> weightedPointMask(4);
1739  std::vector<double> values(8), refValues(8);
1740 
1741 
1742  IntAccessorT idxAcc(mIdxTree);
1743 
1744  // reference data accessors
1745  boost::scoped_ptr<Int16CAccessorT> refSignAcc;
1746  if (mRefSignTreePt) refSignAcc.reset(new Int16CAccessorT(*mRefSignTreePt));
1747 
1748  boost::scoped_ptr<IntCAccessorT> refIdxAcc;
1749  if (mRefIdxTreePt) refIdxAcc.reset(new IntCAccessorT(*mRefIdxTreePt));
1750 
1751  boost::scoped_ptr<AccessorT> refDistAcc;
1752  if (mRefDistTreePt) refDistAcc.reset(new AccessorT(*mRefDistTreePt));
1753 
1754 
1755  for (size_t n = range.begin(); n != range.end(); ++n) {
1756 
1757  coord = mSignLeafs.leaf(n).origin();
1758 
1759  const typename TreeT::LeafNodeType *leafPt = mDistAcc.probeConstLeaf(coord);
1760  typename IntTreeT::LeafNodeType *idxLeafPt = idxAcc.probeLeaf(coord);
1761 
1762 
1763  // reference data leafs
1764  const typename Int16TreeT::LeafNodeType *refSignLeafPt = NULL;
1765  if (refSignAcc) refSignLeafPt = refSignAcc->probeConstLeaf(coord);
1766 
1767  const typename IntTreeT::LeafNodeType *refIdxLeafPt = NULL;
1768  if (refIdxAcc) refIdxLeafPt = refIdxAcc->probeConstLeaf(coord);
1769 
1770  const typename TreeT::LeafNodeType *refDistLeafPt = NULL;
1771  if (refDistAcc) refDistLeafPt = refDistAcc->probeConstLeaf(coord);
1772 
1773 
1774  // generate cell points
1775  size_t ptnIdx = mIndices[n];
1776  coord.offset(TreeT::LeafNodeType::DIM - 1);
1777 
1778 
1779 
1780  for (iter = idxLeafPt->beginValueOn(); iter; ++iter) {
1781 
1782  if(iter.getValue() != 0) continue;
1783 
1784  iter.setValue(ptnIdx);
1785  iter.setValueOff();
1786  offset = iter.pos();
1787  ijk = iter.getCoord();
1788 
1789  const bool inclusiveCell = ijk[0] < coord[0] && ijk[1] < coord[1] && ijk[2] < coord[2];
1790 
1791  const Int16& flags = mSignLeafs.leaf(n).getValue(offset);
1792  signs = SIGNS & flags;
1793  refSigns = 0;
1794 
1795  if ((flags & SEAM) && refSignLeafPt && refIdxLeafPt) {
1796  if (refSignLeafPt->isValueOn(offset)) {
1797  refSigns = (SIGNS & refSignLeafPt->getValue(offset));
1798  }
1799  }
1800 
1801 
1802  if (inclusiveCell) collectCornerValues(*leafPt, offset, values);
1803  else collectCornerValues(mDistAcc, ijk, values);
1804 
1805 
1806  points.clear();
1807  weightedPointMask.clear();
1808 
1809  if (refSigns == 0) {
1810  computeCellPoints(points, values, signs, mIsovalue);
1811  } else {
1812 
1813  if (inclusiveCell) collectCornerValues(*refDistLeafPt, offset, refValues);
1814  else collectCornerValues(*refDistAcc, ijk, refValues);
1815 
1816  computeCellPoints(points, weightedPointMask, values, refValues, signs, refSigns,
1817  mIsovalue, refIdxLeafPt->getValue(offset), *mSeamPointsPt);
1818  }
1819 
1820 
1821  for (size_t i = 0, I = points.size(); i < I; ++i) {
1822 
1823  // offset by cell-origin
1824  points[i][0] += double(ijk[0]);
1825  points[i][1] += double(ijk[1]);
1826  points[i][2] += double(ijk[2]);
1827 
1828 
1829  points[i] = mTransform.indexToWorld(points[i]);
1830 
1831  mPoints[ptnIdx][0] = float(points[i][0]);
1832  mPoints[ptnIdx][1] = float(points[i][1]);
1833  mPoints[ptnIdx][2] = float(points[i][2]);
1834 
1835  if (mSeamPointMaskPt && !weightedPointMask.empty() && weightedPointMask[i]) {
1836  (*mSeamPointMaskPt)[ptnIdx] = 1;
1837  }
1838 
1839  ++ptnIdx;
1840  }
1841  }
1842 
1843  // generate collapsed region points
1844  int onVoxelCount = int(idxLeafPt->onVoxelCount());
1845  while (onVoxelCount > 0) {
1846 
1847  iter = idxLeafPt->beginValueOn();
1848  int regionId = iter.getValue(), count = 0;
1849 
1850  Vec3d avg(0.0), point;
1851 
1852  for (; iter; ++iter) {
1853  if (iter.getValue() != regionId) continue;
1854 
1855  iter.setValue(ptnIdx);
1856  iter.setValueOff();
1857  --onVoxelCount;
1858 
1859  ijk = iter.getCoord();
1860  offset = iter.pos();
1861 
1862  signs = (SIGNS & mSignLeafs.leaf(n).getValue(offset));
1863 
1864  if (ijk[0] < coord[0] && ijk[1] < coord[1] && ijk[2] < coord[2]) {
1865  collectCornerValues(*leafPt, offset, values);
1866  } else {
1867  collectCornerValues(mDistAcc, ijk, values);
1868  }
1869 
1870  points.clear();
1871  computeCellPoints(points, values, signs, mIsovalue);
1872 
1873  avg[0] += double(ijk[0]) + points[0][0];
1874  avg[1] += double(ijk[1]) + points[0][1];
1875  avg[2] += double(ijk[2]) + points[0][2];
1876 
1877  ++count;
1878  }
1879 
1880 
1881  if (count > 1) {
1882  double w = 1.0 / double(count);
1883  avg[0] *= w;
1884  avg[1] *= w;
1885  avg[2] *= w;
1886  }
1887 
1888  avg = mTransform.indexToWorld(avg);
1889 
1890  mPoints[ptnIdx][0] = float(avg[0]);
1891  mPoints[ptnIdx][1] = float(avg[1]);
1892  mPoints[ptnIdx][2] = float(avg[2]);
1893 
1894  ++ptnIdx;
1895  }
1896  }
1897 }
1898 
1899 
1901 
1902 
1903 template<typename TreeT>
1905 {
1906 public:
1908 
1909  typedef typename TreeT::template ValueConverter<int>::Type IntTreeT;
1911 
1912  typedef typename TreeT::template ValueConverter<Int16>::Type Int16TreeT;
1914 
1915  typedef boost::scoped_array<uint32_t> QuantizedPointList;
1916 
1918 
1919  SeamWeights(const TreeT& distTree, const Int16TreeT& refSignTree,
1920  IntTreeT& refIdxTree, QuantizedPointList& points, double iso);
1921 
1922  template <typename LeafNodeType>
1923  void operator()(LeafNodeType &signLeaf, size_t leafIndex) const;
1924 
1925 private:
1926  AccessorT mDistAcc;
1927  Int16AccessorT mRefSignAcc;
1928  IntAccessorT mRefIdxAcc;
1929 
1930  QuantizedPointList& mPoints;
1931  const double mIsovalue;
1932 };
1933 
1934 
1935 template<typename TreeT>
1936 SeamWeights<TreeT>::SeamWeights(const TreeT& distTree, const Int16TreeT& refSignTree,
1937  IntTreeT& refIdxTree, QuantizedPointList& points, double iso)
1938  : mDistAcc(distTree)
1939  , mRefSignAcc(refSignTree)
1940  , mRefIdxAcc(refIdxTree)
1941  , mPoints(points)
1942  , mIsovalue(iso)
1943 {
1944 }
1945 
1946 
1947 template<typename TreeT>
1948 template <typename LeafNodeType>
1949 void
1950 SeamWeights<TreeT>::operator()(LeafNodeType &signLeaf, size_t /*leafIndex*/) const
1951 {
1952  Coord coord = signLeaf.origin();
1953  const typename Int16TreeT::LeafNodeType *refSignLeafPt = mRefSignAcc.probeConstLeaf(coord);
1954 
1955  if (!refSignLeafPt) return;
1956 
1957  const typename TreeT::LeafNodeType *distLeafPt = mDistAcc.probeConstLeaf(coord);
1958  const typename IntTreeT::LeafNodeType *refIdxLeafPt = mRefIdxAcc.probeConstLeaf(coord);
1959 
1960  std::vector<double> values(8);
1961  unsigned char lhsSigns, rhsSigns;
1962  Vec3d point;
1963  Index offset;
1964 
1965  Coord ijk;
1966  coord.offset(TreeT::LeafNodeType::DIM - 1);
1967 
1968  typename LeafNodeType::ValueOnCIter iter = signLeaf.cbeginValueOn();
1969  for (; iter; ++iter) {
1970 
1971  offset = iter.pos();
1972  ijk = iter.getCoord();
1973 
1974  const bool inclusiveCell = ijk[0] < coord[0] && ijk[1] < coord[1] && ijk[2] < coord[2];
1975 
1976  if ((iter.getValue() & SEAM) && refSignLeafPt->isValueOn(offset)) {
1977 
1978  lhsSigns = SIGNS & iter.getValue();
1979  rhsSigns = SIGNS & refSignLeafPt->getValue(offset);
1980 
1981 
1982  if (inclusiveCell) {
1983  collectCornerValues(*distLeafPt, offset, values);
1984  } else {
1985  collectCornerValues(mDistAcc, ijk, values);
1986  }
1987 
1988 
1989  for (size_t n = 1, N = sEdgeGroupTable[lhsSigns][0] + 1; n < N; ++n) {
1990 
1991  int id = matchEdgeGroup(n, lhsSigns, rhsSigns);
1992 
1993  if (id != -1) {
1994 
1995  uint32_t& data = mPoints[refIdxLeafPt->getValue(offset) + (id - 1)];
1996 
1997  if (!(data & MASK_DIRTY_BIT)) {
1998 
1999  int smaples = computeMaskedPoint(point, values, lhsSigns, rhsSigns, n, mIsovalue);
2000 
2001  if (smaples > 0) data = packPoint(point);
2002  else data = MASK_INVALID_BIT;
2003 
2004  data |= MASK_DIRTY_BIT;
2005  }
2006  }
2007  }
2008  }
2009  }
2010 }
2011 
2012 
2014 
2015 
2016 template <typename TreeT, typename LeafManagerT>
2018 {
2019 public:
2020  typedef typename TreeT::ValueType ValueT;
2022 
2023  typedef typename TreeT::template ValueConverter<int>::Type IntTreeT;
2025 
2026  typedef typename TreeT::template ValueConverter<bool>::Type BoolTreeT;
2027 
2028  typedef typename LeafManagerT::TreeType::template ValueConverter<Int16>::Type Int16TreeT;
2030 
2031  typedef typename TreeT::template ValueConverter<float>::Type FloatTreeT;
2033 
2034 
2036 
2037  MergeVoxelRegions(const LeafManagerT& signLeafs, const Int16TreeT& signTree,
2038  const TreeT& distTree, IntTreeT& idxTree, ValueT iso, ValueT adaptivity);
2039 
2040  void run(bool threaded = true);
2041 
2042  void setSpatialAdaptivity(
2043  const math::Transform& distGridXForm, const FloatGridT& adaptivityField);
2044 
2045  void setAdaptivityMask(const BoolTreeT* mask);
2046 
2047  void setRefData(const Int16TreeT* signTree, ValueT adaptivity);
2048 
2050 
2051 
2052  void operator()(const tbb::blocked_range<size_t>&) const;
2053 
2054 private:
2055 
2056  const LeafManagerT& mSignLeafs;
2057 
2058  const Int16TreeT& mSignTree;
2059  Int16AccessorT mSignAcc;
2060 
2061  const TreeT& mDistTree;
2062  AccessorT mDistAcc;
2063 
2064  IntTreeT& mIdxTree;
2065  ValueT mIsovalue, mSurfaceAdaptivity, mInternalAdaptivity;
2066 
2067  const math::Transform* mTransform;
2068  const FloatGridT* mAdaptivityGrid;
2069  const BoolTreeT* mMask;
2070 
2071  const Int16TreeT* mRefSignTree;
2072 };
2073 
2074 template <typename TreeT, typename LeafManagerT>
2076  const LeafManagerT& signLeafs, const Int16TreeT& signTree,
2077  const TreeT& distTree, IntTreeT& idxTree, ValueT iso, ValueT adaptivity)
2078  : mSignLeafs(signLeafs)
2079  , mSignTree(signTree)
2080  , mSignAcc(mSignTree)
2081  , mDistTree(distTree)
2082  , mDistAcc(mDistTree)
2083  , mIdxTree(idxTree)
2084  , mIsovalue(iso)
2085  , mSurfaceAdaptivity(adaptivity)
2086  , mInternalAdaptivity(adaptivity)
2087  , mTransform(NULL)
2088  , mAdaptivityGrid(NULL)
2089  , mMask(NULL)
2090  , mRefSignTree(NULL)
2091 {
2092 }
2093 
2094 
2095 template <typename TreeT, typename LeafManagerT>
2096 void
2098 {
2099  if (threaded) tbb::parallel_for(mSignLeafs.getRange(), *this);
2100  else (*this)(mSignLeafs.getRange());
2101 }
2102 
2103 
2104 template <typename TreeT, typename LeafManagerT>
2105 void
2107  const math::Transform& distGridXForm, const FloatGridT& adaptivityField)
2108 {
2109  mTransform = &distGridXForm;
2110  mAdaptivityGrid = &adaptivityField;
2111 }
2112 
2113 
2114 template <typename TreeT, typename LeafManagerT>
2115 void
2117 {
2118  mMask = mask;
2119 }
2120 
2121 template <typename TreeT, typename LeafManagerT>
2122 void
2124 {
2125  mRefSignTree = signTree;
2126  mInternalAdaptivity = adaptivity;
2127 }
2128 
2129 
2130 template <typename TreeT, typename LeafManagerT>
2131 void
2132 MergeVoxelRegions<TreeT, LeafManagerT>::operator()(const tbb::blocked_range<size_t>& range) const
2133 {
2134  typedef math::Vec3<ValueT> Vec3T;
2135 
2136  typedef typename TreeT::LeafNodeType LeafT;
2137  typedef typename IntTreeT::LeafNodeType IntLeafT;
2138  typedef typename BoolTreeT::LeafNodeType BoolLeafT;
2139  typedef typename LeafT::template ValueConverter<Vec3T>::Type Vec3LeafT;
2140 
2141  const int LeafDim = LeafT::DIM;
2142 
2143  IntAccessorT idxAcc(mIdxTree);
2144 
2145  typename LeafManagerT::TreeType::LeafNodeType::ValueOnCIter iter;
2146 
2147  typedef typename tree::ValueAccessor<const FloatTreeT> FloatTreeCAccessorT;
2148  boost::scoped_ptr<FloatTreeCAccessorT> adaptivityAcc;
2149  if (mAdaptivityGrid) {
2150  adaptivityAcc.reset(new FloatTreeCAccessorT(mAdaptivityGrid->tree()));
2151  }
2152 
2153  typedef typename tree::ValueAccessor<const Int16TreeT> Int16TreeCAccessorT;
2154  boost::scoped_ptr<Int16TreeCAccessorT> refAcc;
2155  if (mRefSignTree) {
2156  refAcc.reset(new Int16TreeCAccessorT(*mRefSignTree));
2157  }
2158 
2159  typedef typename tree::ValueAccessor<const BoolTreeT> BoolTreeCAccessorT;
2160  boost::scoped_ptr<BoolTreeCAccessorT> maskAcc;
2161  if (mMask) {
2162  maskAcc.reset(new BoolTreeCAccessorT(*mMask));
2163  }
2164 
2165  // Allocate reusable leaf buffers
2166  BoolLeafT mask;
2167  Vec3LeafT gradientBuffer;
2168  Coord ijk, nijk, coord, end;
2169 
2170  for (size_t n = range.begin(); n != range.end(); ++n) {
2171 
2172  const Coord& origin = mSignLeafs.leaf(n).origin();
2173 
2174  ValueT adaptivity = mSurfaceAdaptivity;
2175 
2176  if (refAcc && refAcc->probeConstLeaf(origin) == NULL) {
2177  adaptivity = mInternalAdaptivity;
2178  }
2179 
2180  IntLeafT& idxLeaf = *idxAcc.probeLeaf(origin);
2181 
2182  end[0] = origin[0] + LeafDim;
2183  end[1] = origin[1] + LeafDim;
2184  end[2] = origin[2] + LeafDim;
2185 
2186  mask.setValuesOff();
2187 
2188  // Mask off seam line adjacent voxels
2189  if (maskAcc) {
2190  const BoolLeafT* maskLeaf = maskAcc->probeConstLeaf(origin);
2191  if (maskLeaf != NULL) {
2192  typename BoolLeafT::ValueOnCIter it;
2193  for (it = maskLeaf->cbeginValueOn(); it; ++it) {
2194  ijk = it.getCoord();
2195  coord[0] = ijk[0] - (ijk[0] % 2);
2196  coord[1] = ijk[1] - (ijk[1] % 2);
2197  coord[2] = ijk[2] - (ijk[2] % 2);
2198  mask.setActiveState(coord, true);
2199  }
2200  }
2201  }
2202 
2203 
2204  LeafT adaptivityLeaf(origin, adaptivity);
2205 
2206  if (mAdaptivityGrid) {
2207  for (Index offset = 0; offset < LeafT::NUM_VALUES; ++offset) {
2208  ijk = adaptivityLeaf.offsetToGlobalCoord(offset);
2209  Vec3d xyz = mAdaptivityGrid->transform().worldToIndex(
2210  mTransform->indexToWorld(ijk));
2211  ValueT tmpA = ValueT(adaptivityAcc->getValue(util::nearestCoord(xyz)));
2212  adaptivityLeaf.setValueOnly(offset, tmpA * adaptivity);
2213  }
2214  }
2215 
2216  // Mask off ambiguous voxels
2217  for (iter = mSignLeafs.leaf(n).cbeginValueOn(); iter; ++iter) {
2218  ijk = iter.getCoord();
2219  coord[0] = ijk[0] - (ijk[0] % 2);
2220  coord[1] = ijk[1] - (ijk[1] % 2);
2221  coord[2] = ijk[2] - (ijk[2] % 2);
2222  if(mask.isValueOn(coord)) continue;
2223 
2224 
2225 
2226  int flags = int(iter.getValue());
2227  unsigned char signs = SIGNS & flags;
2228  if ((flags & SEAM) || !sAdaptable[signs] || sEdgeGroupTable[signs][0] > 1) {
2229  mask.setActiveState(coord, true);
2230  continue;
2231  }
2232 
2233  for (int i = 0; i < 26; ++i) {
2234  nijk = ijk + util::COORD_OFFSETS[i];
2235  signs = SIGNS & mSignAcc.getValue(nijk);
2236  if (!sAdaptable[signs] || sEdgeGroupTable[signs][0] > 1) {
2237  mask.setActiveState(coord, true);
2238  break;
2239  }
2240  }
2241  }
2242 
2243  int dim = 2;
2244  // Mask off topologically ambiguous 2x2x2 voxel sub-blocks
2245  for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
2246  for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
2247  for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
2248  if (isNonManifold(mDistAcc, ijk, mIsovalue, dim)) {
2249  mask.setActiveState(ijk, true);
2250  }
2251  }
2252  }
2253  }
2254 
2255  // Compute the gradient for the remaining voxels
2256  gradientBuffer.setValuesOff();
2257  for (iter = mSignLeafs.leaf(n).cbeginValueOn(); iter; ++iter) {
2258 
2259  ijk = iter.getCoord();
2260  coord[0] = ijk[0] - (ijk[0] % dim);
2261  coord[1] = ijk[1] - (ijk[1] % dim);
2262  coord[2] = ijk[2] - (ijk[2] % dim);
2263  if(mask.isValueOn(coord)) continue;
2264 
2265  Vec3T norm(math::ISGradient<math::CD_2ND>::result(mDistAcc, ijk));
2266  // Normalize (Vec3's normalize uses isApproxEqual, which uses abs and does more work)
2267  ValueT length = norm.length();
2268  if (length > ValueT(1.0e-7)) {
2269  norm *= ValueT(1.0) / length;
2270  }
2271  gradientBuffer.setValue(ijk, norm);
2272  }
2273 
2274  int regionId = 1, next_dim = dim << 1;
2275 
2276  // Process the first adaptivity level.
2277  for (ijk[0] = 0; ijk[0] < LeafDim; ijk[0] += dim) {
2278  coord[0] = ijk[0] - (ijk[0] % next_dim);
2279  for (ijk[1] = 0; ijk[1] < LeafDim; ijk[1] += dim) {
2280  coord[1] = ijk[1] - (ijk[1] % next_dim);
2281  for (ijk[2] = 0; ijk[2] < LeafDim; ijk[2] += dim) {
2282  coord[2] = ijk[2] - (ijk[2] % next_dim);
2283  adaptivity = adaptivityLeaf.getValue(ijk);
2284  if (mask.isValueOn(ijk) || !isMergable(gradientBuffer, ijk, dim, adaptivity)) {
2285  mask.setActiveState(coord, true);
2286  continue;
2287  }
2288  mergeVoxels(idxLeaf, ijk, dim, regionId++);
2289  }
2290  }
2291  }
2292 
2293 
2294  // Process remaining adaptivity levels
2295  for (dim = 4; dim < LeafDim; dim = dim << 1) {
2296  next_dim = dim << 1;
2297  coord[0] = ijk[0] - (ijk[0] % next_dim);
2298  for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
2299  coord[1] = ijk[1] - (ijk[1] % next_dim);
2300  for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
2301  coord[2] = ijk[2] - (ijk[2] % next_dim);
2302  for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
2303  adaptivity = adaptivityLeaf.getValue(ijk);
2304  if (mask.isValueOn(ijk) || isNonManifold(mDistAcc, ijk, mIsovalue, dim) ||
2305  !isMergable(gradientBuffer, ijk, dim, adaptivity))
2306  {
2307  mask.setActiveState(coord, true);
2308  continue;
2309  }
2310  mergeVoxels(idxLeaf, ijk, dim, regionId++);
2311  }
2312  }
2313  }
2314  }
2315 
2316  adaptivity = adaptivityLeaf.getValue(origin);
2317  if (!(mask.isValueOn(origin) || isNonManifold(mDistAcc, origin, mIsovalue, LeafDim))
2318  && isMergable(gradientBuffer, origin, LeafDim, adaptivity))
2319  {
2320  mergeVoxels(idxLeaf, origin, LeafDim, regionId++);
2321  }
2322  }
2323 }
2324 
2325 
2327 
2328 
2329 // Constructs qudas
2331 {
2332  UniformPrimBuilder(): mIdx(0), mPolygonPool(NULL) {}
2333 
2334  void init(const size_t upperBound, PolygonPool& quadPool)
2335  {
2336  mPolygonPool = &quadPool;
2337  mPolygonPool->resetQuads(upperBound);
2338  mIdx = 0;
2339  }
2340 
2341  void addPrim(const Vec4I& verts, bool reverse, char flags = 0)
2342  {
2343  if (!reverse) {
2344  mPolygonPool->quad(mIdx) = verts;
2345  } else {
2346  Vec4I& quad = mPolygonPool->quad(mIdx);
2347  quad[0] = verts[3];
2348  quad[1] = verts[2];
2349  quad[2] = verts[1];
2350  quad[3] = verts[0];
2351  }
2352  mPolygonPool->quadFlags(mIdx) = flags;
2353  ++mIdx;
2354  }
2355 
2356  void done()
2357  {
2358  mPolygonPool->trimQuads(mIdx);
2359  }
2360 
2361 private:
2362  size_t mIdx;
2363  PolygonPool* mPolygonPool;
2364 };
2365 
2366 
2367 // Constructs qudas and triangles
2369 {
2370  AdaptivePrimBuilder() : mQuadIdx(0), mTriangleIdx(0), mPolygonPool(NULL) {}
2371 
2372  void init(const size_t upperBound, PolygonPool& polygonPool)
2373  {
2374  mPolygonPool = &polygonPool;
2375  mPolygonPool->resetQuads(upperBound);
2376  mPolygonPool->resetTriangles(upperBound);
2377 
2378  mQuadIdx = 0;
2379  mTriangleIdx = 0;
2380  }
2381 
2382  void addPrim(const Vec4I& verts, bool reverse, char flags = 0)
2383  {
2384  if (verts[0] != verts[1] && verts[0] != verts[2] && verts[0] != verts[3]
2385  && verts[1] != verts[2] && verts[1] != verts[3] && verts[2] != verts[3]) {
2386  mPolygonPool->quadFlags(mQuadIdx) = flags;
2387  addQuad(verts, reverse);
2388  } else if (
2389  verts[0] == verts[3] &&
2390  verts[1] != verts[2] &&
2391  verts[1] != verts[0] &&
2392  verts[2] != verts[0]) {
2393  mPolygonPool->triangleFlags(mTriangleIdx) = flags;
2394  addTriangle(verts[0], verts[1], verts[2], reverse);
2395  } else if (
2396  verts[1] == verts[2] &&
2397  verts[0] != verts[3] &&
2398  verts[0] != verts[1] &&
2399  verts[3] != verts[1]) {
2400  mPolygonPool->triangleFlags(mTriangleIdx) = flags;
2401  addTriangle(verts[0], verts[1], verts[3], reverse);
2402  } else if (
2403  verts[0] == verts[1] &&
2404  verts[2] != verts[3] &&
2405  verts[2] != verts[0] &&
2406  verts[3] != verts[0]) {
2407  mPolygonPool->triangleFlags(mTriangleIdx) = flags;
2408  addTriangle(verts[0], verts[2], verts[3], reverse);
2409  } else if (
2410  verts[2] == verts[3] &&
2411  verts[0] != verts[1] &&
2412  verts[0] != verts[2] &&
2413  verts[1] != verts[2]) {
2414  mPolygonPool->triangleFlags(mTriangleIdx) = flags;
2415  addTriangle(verts[0], verts[1], verts[2], reverse);
2416  }
2417  }
2418 
2419 
2420  void done()
2421  {
2422  mPolygonPool->trimQuads(mQuadIdx, /*reallocate=*/true);
2423  mPolygonPool->trimTrinagles(mTriangleIdx, /*reallocate=*/true);
2424  }
2425 
2426 private:
2427 
2428  void addQuad(const Vec4I& verts, bool reverse)
2429  {
2430  if (!reverse) {
2431  mPolygonPool->quad(mQuadIdx) = verts;
2432  } else {
2433  Vec4I& quad = mPolygonPool->quad(mQuadIdx);
2434  quad[0] = verts[3];
2435  quad[1] = verts[2];
2436  quad[2] = verts[1];
2437  quad[3] = verts[0];
2438  }
2439  ++mQuadIdx;
2440  }
2441 
2442  void addTriangle(unsigned v0, unsigned v1, unsigned v2, bool reverse)
2443  {
2444  Vec3I& prim = mPolygonPool->triangle(mTriangleIdx);
2445 
2446  prim[1] = v1;
2447 
2448  if (!reverse) {
2449  prim[0] = v0;
2450  prim[2] = v2;
2451  } else {
2452  prim[0] = v2;
2453  prim[2] = v0;
2454  }
2455  ++mTriangleIdx;
2456  }
2457 
2458  size_t mQuadIdx, mTriangleIdx;
2459  PolygonPool *mPolygonPool;
2460 };
2461 
2462 
2463 template<typename SignAccT, typename IdxAccT, typename PrimBuilder>
2464 inline void
2465 constructPolygons(Int16 flags, Int16 refFlags, const Vec4i& offsets, const Coord& ijk,
2466  const SignAccT& signAcc, const IdxAccT& idxAcc, PrimBuilder& mesher, Index32 pointListSize)
2467 {
2468  const Index32 v0 = idxAcc.getValue(ijk);
2469  if (v0 == util::INVALID_IDX) return;
2470 
2471  char tag[2];
2472  tag[0] = (flags & SEAM) ? POLYFLAG_FRACTURE_SEAM : 0;
2473  tag[1] = tag[0] | char(POLYFLAG_EXTERIOR);
2474 
2475  const bool isInside = flags & INSIDE;
2476 
2477  Coord coord;
2478  openvdb::Vec4I quad;
2479  unsigned char cell;
2480  Index32 tmpIdx = 0;
2481 
2482  if (flags & XEDGE) {
2483 
2484  quad[0] = v0 + offsets[0];
2485 
2486  // i, j-1, k
2487  coord[0] = ijk[0];
2488  coord[1] = ijk[1] - 1;
2489  coord[2] = ijk[2];
2490 
2491  quad[1] = idxAcc.getValue(coord);
2492  cell = SIGNS & signAcc.getValue(coord);
2493  if (sEdgeGroupTable[cell][0] > 1) {
2494  tmpIdx = quad[1] + Index32(sEdgeGroupTable[cell][5] - 1);
2495  if (tmpIdx < pointListSize) quad[1] = tmpIdx;
2496  }
2497 
2498  // i, j-1, k-1
2499  coord[2] -= 1;
2500 
2501  quad[2] = idxAcc.getValue(coord);
2502  cell = SIGNS & signAcc.getValue(coord);
2503  if (sEdgeGroupTable[cell][0] > 1) {
2504  tmpIdx = quad[2] + Index32(sEdgeGroupTable[cell][7] - 1);
2505  if (tmpIdx < pointListSize) quad[2] = tmpIdx;
2506  }
2507 
2508  // i, j, k-1
2509  coord[1] = ijk[1];
2510 
2511  quad[3] = idxAcc.getValue(coord);
2512  cell = SIGNS & signAcc.getValue(coord);
2513  if (sEdgeGroupTable[cell][0] > 1) {
2514  tmpIdx = quad[3] + Index32(sEdgeGroupTable[cell][3] - 1);
2515  if (tmpIdx < pointListSize) quad[3] = tmpIdx;
2516  }
2517 
2518  if (quad[1] != util::INVALID_IDX &&
2519  quad[2] != util::INVALID_IDX && quad[3] != util::INVALID_IDX) {
2520  mesher.addPrim(quad, isInside, tag[bool(refFlags & XEDGE)]);
2521  }
2522  }
2523 
2524 
2525  if (flags & YEDGE) {
2526 
2527  quad[0] = v0 + offsets[1];
2528 
2529  // i, j, k-1
2530  coord[0] = ijk[0];
2531  coord[1] = ijk[1];
2532  coord[2] = ijk[2] - 1;
2533 
2534  quad[1] = idxAcc.getValue(coord);
2535  cell = SIGNS & signAcc.getValue(coord);
2536  if (sEdgeGroupTable[cell][0] > 1) {
2537  tmpIdx = quad[1] + Index32(sEdgeGroupTable[cell][12] - 1);
2538  if (tmpIdx < pointListSize) quad[1] = tmpIdx;
2539  }
2540 
2541  // i-1, j, k-1
2542  coord[0] -= 1;
2543 
2544  quad[2] = idxAcc.getValue(coord);
2545  cell = SIGNS & signAcc.getValue(coord);
2546  if (sEdgeGroupTable[cell][0] > 1) {
2547  tmpIdx = quad[2] + Index32(sEdgeGroupTable[cell][11] - 1);
2548  if (tmpIdx < pointListSize) quad[2] = tmpIdx;
2549  }
2550 
2551  // i-1, j, k
2552  coord[2] = ijk[2];
2553 
2554  quad[3] = idxAcc.getValue(coord);
2555  cell = SIGNS & signAcc.getValue(coord);
2556  if (sEdgeGroupTable[cell][0] > 1) {
2557  tmpIdx = quad[3] + Index32(sEdgeGroupTable[cell][10] - 1);
2558  if (tmpIdx < pointListSize) quad[3] = tmpIdx;
2559  }
2560 
2561  if (quad[1] != util::INVALID_IDX &&
2562  quad[2] != util::INVALID_IDX && quad[3] != util::INVALID_IDX) {
2563  mesher.addPrim(quad, isInside, tag[bool(refFlags & YEDGE)]);
2564  }
2565  }
2566 
2567  if (flags & ZEDGE) {
2568 
2569  quad[0] = v0 + offsets[2];
2570 
2571  // i, j-1, k
2572  coord[0] = ijk[0];
2573  coord[1] = ijk[1] - 1;
2574  coord[2] = ijk[2];
2575 
2576  quad[1] = idxAcc.getValue(coord);
2577  cell = SIGNS & signAcc.getValue(coord);
2578  if (sEdgeGroupTable[cell][0] > 1) {
2579  tmpIdx = quad[1] + Index32(sEdgeGroupTable[cell][8] - 1);
2580  if (tmpIdx < pointListSize) quad[1] = tmpIdx;
2581  }
2582 
2583  // i-1, j-1, k
2584  coord[0] -= 1;
2585 
2586  quad[2] = idxAcc.getValue(coord);
2587  cell = SIGNS & signAcc.getValue(coord);
2588  if (sEdgeGroupTable[cell][0] > 1) {
2589  tmpIdx = quad[2] + Index32(sEdgeGroupTable[cell][6] - 1);
2590  if (tmpIdx < pointListSize) quad[2] = tmpIdx;
2591  }
2592 
2593  // i-1, j, k
2594  coord[1] = ijk[1];
2595 
2596  quad[3] = idxAcc.getValue(coord);
2597  cell = SIGNS & signAcc.getValue(coord);
2598  if (sEdgeGroupTable[cell][0] > 1) {
2599  tmpIdx = quad[3] + Index32(sEdgeGroupTable[cell][2] - 1);
2600  if (tmpIdx < pointListSize) quad[3] = tmpIdx;
2601  }
2602 
2603  if (quad[1] != util::INVALID_IDX &&
2604  quad[2] != util::INVALID_IDX && quad[3] != util::INVALID_IDX) {
2605  mesher.addPrim(quad, !isInside, tag[bool(refFlags & ZEDGE)]);
2606  }
2607  }
2608 }
2609 
2610 
2612 
2613 
2614 template<typename LeafManagerT, typename PrimBuilder>
2616 {
2617 public:
2618  typedef typename LeafManagerT::TreeType::template ValueConverter<int>::Type IntTreeT;
2619  typedef typename LeafManagerT::TreeType::template ValueConverter<Int16>::Type Int16TreeT;
2620 
2623 
2625 
2626 
2627  GenPolygons(const LeafManagerT& signLeafs, const Int16TreeT& signTree,
2628  const IntTreeT& idxTree, PolygonPoolList& polygons, Index32 pointListSize);
2629 
2630  void run(bool threaded = true);
2631 
2632 
2633  void setRefSignTree(const Int16TreeT *r) { mRefSignTree = r; }
2634 
2636 
2637 
2638  void operator()(const tbb::blocked_range<size_t>&) const;
2639 
2640 private:
2641  const LeafManagerT& mSignLeafs;
2642  const Int16TreeT& mSignTree;
2643  const IntTreeT& mIdxTree;
2644  const PolygonPoolList& mPolygonPoolList;
2645  const Index32 mPointListSize;
2646 
2647  const Int16TreeT *mRefSignTree;
2648  };
2649 
2650 
2651 template<typename LeafManagerT, typename PrimBuilder>
2653  const Int16TreeT& signTree, const IntTreeT& idxTree, PolygonPoolList& polygons,
2654  Index32 pointListSize)
2655  : mSignLeafs(signLeafs)
2656  , mSignTree(signTree)
2657  , mIdxTree(idxTree)
2658  , mPolygonPoolList(polygons)
2659  , mPointListSize(pointListSize)
2660  , mRefSignTree(NULL)
2661 {
2662 }
2663 
2664 template<typename LeafManagerT, typename PrimBuilder>
2665 void
2667 {
2668  if (threaded) tbb::parallel_for(mSignLeafs.getRange(), *this);
2669  else (*this)(mSignLeafs.getRange());
2670 }
2671 
2672 template<typename LeafManagerT, typename PrimBuilder>
2673 void
2675  const tbb::blocked_range<size_t>& range) const
2676 {
2677  typename LeafManagerT::TreeType::LeafNodeType::ValueOnCIter iter;
2678  IntAccessorT idxAcc(mIdxTree);
2679  Int16AccessorT signAcc(mSignTree);
2680 
2681 
2682  PrimBuilder mesher;
2683  size_t edgeCount;
2684  Coord ijk, origin;
2685 
2686 
2687  // reference data
2688  boost::scoped_ptr<Int16AccessorT> refSignAcc;
2689  if (mRefSignTree) refSignAcc.reset(new Int16AccessorT(*mRefSignTree));
2690 
2691 
2692  for (size_t n = range.begin(); n != range.end(); ++n) {
2693 
2694  origin = mSignLeafs.leaf(n).origin();
2695 
2696  // Get an upper bound on the number of primitives.
2697  edgeCount = 0;
2698  iter = mSignLeafs.leaf(n).cbeginValueOn();
2699  for (; iter; ++iter) {
2700  if (iter.getValue() & XEDGE) ++edgeCount;
2701  if (iter.getValue() & YEDGE) ++edgeCount;
2702  if (iter.getValue() & ZEDGE) ++edgeCount;
2703  }
2704 
2705  if(edgeCount == 0) continue;
2706 
2707  mesher.init(edgeCount, mPolygonPoolList[n]);
2708 
2709  const typename Int16TreeT::LeafNodeType *signleafPt = signAcc.probeConstLeaf(origin);
2710  const typename IntTreeT::LeafNodeType *idxLeafPt = idxAcc.probeConstLeaf(origin);
2711 
2712  if (!signleafPt || !idxLeafPt) continue;
2713 
2714 
2715  const typename Int16TreeT::LeafNodeType *refSignLeafPt = NULL;
2716  if (refSignAcc) refSignLeafPt = refSignAcc->probeConstLeaf(origin);
2717 
2718  Vec4i offsets;
2719 
2720  iter = mSignLeafs.leaf(n).cbeginValueOn();
2721  for (; iter; ++iter) {
2722  ijk = iter.getCoord();
2723 
2724  Int16 flags = iter.getValue();
2725 
2726  if (!(flags & 0xE00)) continue;
2727 
2728  Int16 refFlags = 0;
2729  if (refSignLeafPt) {
2730  refFlags = refSignLeafPt->getValue(iter.pos());
2731  }
2732 
2733  offsets[0] = 0;
2734  offsets[1] = 0;
2735  offsets[2] = 0;
2736 
2737  const unsigned char cell = (SIGNS & flags);
2738 
2739  if (sEdgeGroupTable[cell][0] > 1) {
2740  offsets[0] = (sEdgeGroupTable[cell][1] - 1);
2741  offsets[1] = (sEdgeGroupTable[cell][9] - 1);
2742  offsets[2] = (sEdgeGroupTable[cell][4] - 1);
2743  }
2744 
2745  if (ijk[0] > origin[0] && ijk[1] > origin[1] && ijk[2] > origin[2]) {
2746  constructPolygons(flags, refFlags, offsets, ijk, *signleafPt, *idxLeafPt, mesher, mPointListSize);
2747  } else {
2748  constructPolygons(flags, refFlags, offsets, ijk, signAcc, idxAcc, mesher, mPointListSize);
2749  }
2750  }
2751 
2752  mesher.done();
2753  }
2754 }
2755 
2756 
2758 
2759 // Masking and mesh partitioning
2760 
2761 struct PartOp
2762 {
2763 
2764  PartOp(size_t leafCount, size_t partitions, size_t activePart)
2765  {
2766  size_t leafSegments = leafCount / partitions;
2767  mStart = leafSegments * activePart;
2768  mEnd = activePart >= (partitions - 1) ? leafCount : mStart + leafSegments;
2769  }
2770 
2771  template <typename LeafNodeType>
2772  void operator()(LeafNodeType &leaf, size_t leafIndex) const
2773  {
2774  if (leafIndex < mStart || leafIndex >= mEnd) leaf.setValuesOff();
2775  }
2776 
2777 private:
2778  size_t mStart, mEnd;
2779 };
2780 
2781 
2783 
2784 
2785 template<typename SrcTreeT>
2786 class PartGen
2787 {
2788 public:
2790  typedef typename SrcTreeT::template ValueConverter<bool>::Type BoolTreeT;
2792 
2794 
2795 
2796  PartGen(const LeafManagerT& leafs, size_t partitions, size_t activePart);
2797 
2798  void run(bool threaded = true);
2799 
2800  BoolTreeT& tree() { return mTree; }
2801 
2802 
2804 
2805  PartGen(PartGen&, tbb::split);
2806  void operator()(const tbb::blocked_range<size_t>&);
2807  void join(PartGen& rhs) { mTree.merge(rhs.mTree); }
2808 
2809 private:
2810  const LeafManagerT& mLeafManager;
2811  BoolTreeT mTree;
2812  size_t mStart, mEnd;
2813 };
2814 
2815 template<typename SrcTreeT>
2816 PartGen<SrcTreeT>::PartGen(const LeafManagerT& leafs, size_t partitions, size_t activePart)
2817  : mLeafManager(leafs)
2818  , mTree(false)
2819  , mStart(0)
2820  , mEnd(0)
2821 {
2822  size_t leafCount = leafs.leafCount();
2823  size_t leafSegments = leafCount / partitions;
2824  mStart = leafSegments * activePart;
2825  mEnd = activePart >= (partitions - 1) ? leafCount : mStart + leafSegments;
2826 }
2827 
2828 template<typename SrcTreeT>
2830  : mLeafManager(rhs.mLeafManager)
2831  , mTree(false)
2832  , mStart(rhs.mStart)
2833  , mEnd(rhs.mEnd)
2834 {
2835 }
2836 
2837 
2838 template<typename SrcTreeT>
2839 void
2841 {
2842  if (threaded) tbb::parallel_reduce(mLeafManager.getRange(), *this);
2843  else (*this)(mLeafManager.getRange());
2844 }
2845 
2846 
2847 template<typename SrcTreeT>
2848 void
2849 PartGen<SrcTreeT>::operator()(const tbb::blocked_range<size_t>& range)
2850 {
2851  Coord ijk;
2852  BoolAccessorT acc(mTree);
2853 
2854  typedef typename BoolTreeT::LeafNodeType BoolLeafT;
2855  typename SrcTreeT::LeafNodeType::ValueOnCIter iter;
2856 
2857  for (size_t n = range.begin(); n != range.end(); ++n) {
2858  if (n < mStart || n >= mEnd) continue;
2859  BoolLeafT* leaf = acc.touchLeaf(mLeafManager.leaf(n).origin());
2860  leaf->topologyUnion(mLeafManager.leaf(n));
2861  }
2862 }
2863 
2864 
2866 
2867 
2868 template<typename TreeT, typename LeafManagerT>
2870 {
2871 public:
2872  typedef typename TreeT::template ValueConverter<bool>::Type BoolTreeT;
2873 
2875 
2876  GenSeamMask(const LeafManagerT& leafs, const TreeT& tree);
2877 
2878  void run(bool threaded = true);
2879 
2880  BoolTreeT& mask() { return mMaskTree; }
2881 
2883 
2884  GenSeamMask(GenSeamMask&, tbb::split);
2885  void operator()(const tbb::blocked_range<size_t>&);
2886  void join(GenSeamMask& rhs) { mMaskTree.merge(rhs.mMaskTree); }
2887 
2888 private:
2889 
2890  const LeafManagerT& mLeafManager;
2891  const TreeT& mTree;
2892 
2893  BoolTreeT mMaskTree;
2894 };
2895 
2896 
2897 template<typename TreeT, typename LeafManagerT>
2898 GenSeamMask<TreeT, LeafManagerT>::GenSeamMask(const LeafManagerT& leafs, const TreeT& tree)
2899  : mLeafManager(leafs)
2900  , mTree(tree)
2901  , mMaskTree(false)
2902 {
2903 }
2904 
2905 
2906 template<typename TreeT, typename LeafManagerT>
2908  : mLeafManager(rhs.mLeafManager)
2909  , mTree(rhs.mTree)
2910  , mMaskTree(false)
2911 {
2912 }
2913 
2914 
2915 template<typename TreeT, typename LeafManagerT>
2916 void
2918 {
2919  if (threaded) tbb::parallel_reduce(mLeafManager.getRange(), *this);
2920  else (*this)(mLeafManager.getRange());
2921 }
2922 
2923 
2924 template<typename TreeT, typename LeafManagerT>
2925 void
2926 GenSeamMask<TreeT, LeafManagerT>::operator()(const tbb::blocked_range<size_t>& range)
2927 {
2928  Coord ijk;
2930  tree::ValueAccessor<BoolTreeT> maskAcc(mMaskTree);
2931 
2932  typename LeafManagerT::TreeType::LeafNodeType::ValueOnCIter it;
2933 
2934  for (size_t n = range.begin(); n != range.end(); ++n) {
2935 
2936  it = mLeafManager.leaf(n).cbeginValueOn();
2937 
2938  for (; it; ++it) {
2939 
2940  ijk = it.getCoord();
2941 
2942  unsigned char rhsSigns = acc.getValue(ijk) & SIGNS;
2943 
2944  if (sEdgeGroupTable[rhsSigns][0] > 0) {
2945  unsigned char lhsSigns = it.getValue() & SIGNS;
2946  if (rhsSigns != lhsSigns) {
2947  maskAcc.setValueOn(ijk);
2948  }
2949  }
2950  }
2951  }
2952 }
2953 
2954 
2956 
2957 
2958 template<typename TreeT>
2960 {
2961 public:
2963 
2964  TagSeamEdges(const TreeT& tree) : mAcc(tree) {}
2965 
2966  template <typename LeafNodeType>
2967  void operator()(LeafNodeType &leaf, size_t/*leafIndex*/) const
2968  {
2969  const typename TreeT::LeafNodeType *maskLeaf =
2970  mAcc.probeConstLeaf(leaf.origin());
2971 
2972  if (!maskLeaf) return;
2973 
2974  typename LeafNodeType::ValueOnIter it = leaf.beginValueOn();
2975 
2976  for (; it; ++it) {
2977 
2978  if (maskLeaf->isValueOn(it.pos())) {
2979  it.setValue(it.getValue() | SEAM);
2980  }
2981  }
2982  }
2983 
2984 private:
2985  AccessorT mAcc;
2986 };
2987 
2988 
2989 
2990 template<typename BoolTreeT>
2992 {
2994 
2995  MaskEdges(const BoolTreeT& valueMask) : mMaskAcc(valueMask) {}
2996 
2997  template <typename LeafNodeType>
2998  void operator()(LeafNodeType &leaf, size_t /*leafIndex*/) const
2999  {
3000  typename LeafNodeType::ValueOnIter it = leaf.beginValueOn();
3001 
3002  const typename BoolTreeT::LeafNodeType * maskLeaf =
3003  mMaskAcc.probeConstLeaf(leaf.origin());
3004 
3005  if (maskLeaf) {
3006  for (; it; ++it) {
3007  if (!maskLeaf->isValueOn(it.pos())) {
3008  it.setValue(0x1FF & it.getValue());
3009  }
3010  }
3011  } else {
3012  for (; it; ++it) {
3013  it.setValue(0x1FF & it.getValue());
3014  }
3015  }
3016  }
3017 
3018 private:
3019  BoolAccessorT mMaskAcc;
3020 };
3021 
3022 
3024 {
3025 public:
3027 
3028  FlagUsedPoints(const PolygonPoolList& polygons, size_t polyListCount,
3029  std::vector<unsigned char>& usedPointMask)
3030  : mPolygons(polygons)
3031  , mPolyListCount(polyListCount)
3032  , mUsedPointMask(usedPointMask)
3033  {
3034  }
3035 
3036  void run(bool threaded = true)
3037  {
3038  if (threaded) {
3039  tbb::parallel_for(tbb::blocked_range<size_t>(0, mPolyListCount), *this);
3040  } else {
3041  (*this)(tbb::blocked_range<size_t>(0, mPolyListCount));
3042  }
3043  }
3044 
3046 
3047  void operator()(const tbb::blocked_range<size_t>& range) const
3048  {
3049  // Concurrent writes to same memory address can occur, but
3050  // all threads are writing the same value and char is atomic.
3051  for (size_t n = range.begin(); n != range.end(); ++n) {
3052  const PolygonPool& polygons = mPolygons[n];
3053  for (size_t i = 0; i < polygons.numQuads(); ++i) {
3054  const Vec4I& quad = polygons.quad(i);
3055  mUsedPointMask[quad[0]] = 1;
3056  mUsedPointMask[quad[1]] = 1;
3057  mUsedPointMask[quad[2]] = 1;
3058  mUsedPointMask[quad[3]] = 1;
3059  }
3060 
3061  for (size_t i = 0; i < polygons.numTriangles(); ++i) {
3062  const Vec3I& triangle = polygons.triangle(i);
3063  mUsedPointMask[triangle[0]] = 1;
3064  mUsedPointMask[triangle[1]] = 1;
3065  mUsedPointMask[triangle[2]] = 1;
3066  }
3067  }
3068  }
3069 
3070 
3071 private:
3072  const PolygonPoolList& mPolygons;
3073  size_t mPolyListCount;
3074  std::vector<unsigned char>& mUsedPointMask;
3075 };
3076 
3078 {
3079 public:
3081 
3083  size_t polyListCount, const std::vector<unsigned>& indexMap)
3084  : mPolygons(polygons)
3085  , mPolyListCount(polyListCount)
3086  , mIndexMap(indexMap)
3087  {
3088  }
3089 
3090  void run(bool threaded = true)
3091  {
3092  if (threaded) {
3093  tbb::parallel_for(tbb::blocked_range<size_t>(0, mPolyListCount), *this);
3094  } else {
3095  (*this)(tbb::blocked_range<size_t>(0, mPolyListCount));
3096  }
3097  }
3098 
3100 
3101  void operator()(const tbb::blocked_range<size_t>& range) const
3102  {
3103  for (size_t n = range.begin(); n != range.end(); ++n) {
3104  PolygonPool& polygons = mPolygons[n];
3105  for (size_t i = 0; i < polygons.numQuads(); ++i) {
3106  Vec4I& quad = polygons.quad(i);
3107  quad[0] = mIndexMap[quad[0]];
3108  quad[1] = mIndexMap[quad[1]];
3109  quad[2] = mIndexMap[quad[2]];
3110  quad[3] = mIndexMap[quad[3]];
3111  }
3112 
3113  for (size_t i = 0; i < polygons.numTriangles(); ++i) {
3114  Vec3I& triangle = polygons.triangle(i);
3115  triangle[0] = mIndexMap[triangle[0]];
3116  triangle[1] = mIndexMap[triangle[1]];
3117  triangle[2] = mIndexMap[triangle[2]];
3118  }
3119  }
3120  }
3121 
3122 
3123 private:
3124  PolygonPoolList& mPolygons;
3125  size_t mPolyListCount;
3126  const std::vector<unsigned>& mIndexMap;
3127 };
3128 
3129 
3131 {
3132 public:
3134 
3137  const PointList& oldPointList,
3138  const std::vector<unsigned>& indexMap,
3139  const std::vector<unsigned char>& usedPointMask)
3140  : mNewPointList(newPointList)
3141  , mOldPointList(oldPointList)
3142  , mIndexMap(indexMap)
3143  , mUsedPointMask(usedPointMask)
3144  {
3145  }
3146 
3147  void run(bool threaded = true)
3148  {
3149  if (threaded) {
3150  tbb::parallel_for(tbb::blocked_range<size_t>(0, mIndexMap.size()), *this);
3151  } else {
3152  (*this)(tbb::blocked_range<size_t>(0, mIndexMap.size()));
3153  }
3154  }
3155 
3157 
3158  void operator()(const tbb::blocked_range<size_t>& range) const
3159  {
3160  for (size_t n = range.begin(); n != range.end(); ++n) {
3161  if (mUsedPointMask[n]) {
3162  const size_t index = mIndexMap[n];
3163  mNewPointList.get()[index] = mOldPointList[n];
3164  }
3165  }
3166  }
3167 
3168 private:
3170  const PointList& mOldPointList;
3171  const std::vector<unsigned>& mIndexMap;
3172  const std::vector<unsigned char>& mUsedPointMask;
3173 };
3174 
3175 
3177 
3178 
3179 template<typename SrcTreeT>
3181 {
3182 public:
3184  typedef typename SrcTreeT::template ValueConverter<bool>::Type BoolTreeT;
3188 
3189 
3191 
3192 
3193  GenTopologyMask(const BoolGridT& mask, const LeafManagerT& srcLeafs,
3194  const math::Transform& srcXForm, bool invertMask);
3195 
3196  void run(bool threaded = true);
3197 
3198  BoolTreeT& tree() { return mTree; }
3199 
3200 
3202 
3203  GenTopologyMask(GenTopologyMask&, tbb::split);
3204 
3205  void operator()(const tbb::blocked_range<size_t>&);
3206 
3207  void join(GenTopologyMask& rhs) { mTree.merge(rhs.mTree); }
3208 
3209 private:
3210 
3211  const BoolGridT& mMask;
3212  const LeafManagerT& mLeafManager;
3213  const math::Transform& mSrcXForm;
3214  bool mInvertMask;
3215  BoolTreeT mTree;
3216 };
3217 
3218 
3219 template<typename SrcTreeT>
3221  const math::Transform& srcXForm, bool invertMask)
3222  : mMask(mask)
3223  , mLeafManager(srcLeafs)
3224  , mSrcXForm(srcXForm)
3225  , mInvertMask(invertMask)
3226  , mTree(false)
3227 {
3228 }
3229 
3230 
3231 template<typename SrcTreeT>
3233  : mMask(rhs.mMask)
3234  , mLeafManager(rhs.mLeafManager)
3235  , mSrcXForm(rhs.mSrcXForm)
3236  , mInvertMask(rhs.mInvertMask)
3237  , mTree(false)
3238 {
3239 }
3240 
3241 
3242 template<typename SrcTreeT>
3243 void
3245 {
3246  if (threaded) {
3247  tbb::parallel_reduce(mLeafManager.getRange(), *this);
3248  } else {
3249  (*this)(mLeafManager.getRange());
3250  }
3251 }
3252 
3253 
3254 template<typename SrcTreeT>
3255 void
3256 GenTopologyMask<SrcTreeT>::operator()(const tbb::blocked_range<size_t>& range)
3257 {
3258  Coord ijk;
3259  Vec3d xyz;
3260  typedef typename BoolTreeT::LeafNodeType BoolLeafT;
3261  const math::Transform& maskXForm = mMask.transform();
3262  tree::ValueAccessor<const BoolTreeT> maskAcc(mMask.tree());
3263  tree::ValueAccessor<BoolTreeT> acc(mTree);
3264 
3265  typename SrcTreeT::LeafNodeType::ValueOnCIter iter;
3266  for (size_t n = range.begin(); n != range.end(); ++n) {
3267 
3268  ijk = mLeafManager.leaf(n).origin();
3269  BoolLeafT* leaf = new BoolLeafT(ijk, false);
3270  bool addLeaf = false;
3271 
3272  if (maskXForm == mSrcXForm) {
3273 
3274  const BoolLeafT* maskLeaf = maskAcc.probeConstLeaf(ijk);
3275 
3276  if (maskLeaf) {
3277 
3278  for (iter = mLeafManager.leaf(n).cbeginValueOn(); iter; ++iter) {
3279  Index pos = iter.pos();
3280  if(maskLeaf->isValueOn(pos) != mInvertMask) {
3281  leaf->setValueOn(pos);
3282  addLeaf = true;
3283  }
3284  }
3285 
3286  } else if (maskAcc.isValueOn(ijk) != mInvertMask) {
3287  leaf->topologyUnion(mLeafManager.leaf(n));
3288  addLeaf = true;
3289  }
3290 
3291  } else {
3292  for (iter = mLeafManager.leaf(n).cbeginValueOn(); iter; ++iter) {
3293  ijk = iter.getCoord();
3294  xyz = maskXForm.worldToIndex(mSrcXForm.indexToWorld(ijk));
3295  if(maskAcc.isValueOn(util::nearestCoord(xyz)) != mInvertMask) {
3296  leaf->setValueOn(iter.pos());
3297  addLeaf = true;
3298  }
3299  }
3300  }
3301 
3302  if (addLeaf) acc.addLeaf(leaf);
3303  else delete leaf;
3304  }
3305 }
3306 
3307 
3309 
3310 
3311 template<typename SrcTreeT>
3313 {
3314 public:
3315  typedef typename SrcTreeT::template ValueConverter<int>::Type IntTreeT;
3316  typedef typename SrcTreeT::template ValueConverter<bool>::Type BoolTreeT;
3318 
3320 
3321  GenBoundaryMask(const LeafManagerT& leafs, const BoolTreeT&, const IntTreeT&);
3322 
3323  void run(bool threaded = true);
3324 
3325  BoolTreeT& tree() { return mTree; }
3326 
3328 
3329  GenBoundaryMask(GenBoundaryMask&, tbb::split);
3330  void operator()(const tbb::blocked_range<size_t>&);
3331  void join(GenBoundaryMask& rhs) { mTree.merge(rhs.mTree); }
3332 
3333 private:
3334  // This typedef is needed for Windows
3335  typedef tree::ValueAccessor<const IntTreeT> IntTreeAccessorT;
3336 
3337  bool neighboringLeaf(const Coord&, const IntTreeAccessorT&) const;
3338 
3339  const LeafManagerT& mLeafManager;
3340  const BoolTreeT& mMaskTree;
3341  const IntTreeT& mIdxTree;
3342  BoolTreeT mTree;
3343  CoordBBox mLeafBBox;
3344 };
3345 
3346 
3347 template<typename SrcTreeT>
3349  const BoolTreeT& maskTree, const IntTreeT& auxTree)
3350  : mLeafManager(leafs)
3351  , mMaskTree(maskTree)
3352  , mIdxTree(auxTree)
3353  , mTree(false)
3354 {
3355  mIdxTree.evalLeafBoundingBox(mLeafBBox);
3356  mLeafBBox.expand(IntTreeT::LeafNodeType::DIM);
3357 }
3358 
3359 
3360 template<typename SrcTreeT>
3362  : mLeafManager(rhs.mLeafManager)
3363  , mMaskTree(rhs.mMaskTree)
3364  , mIdxTree(rhs.mIdxTree)
3365  , mTree(false)
3366  , mLeafBBox(rhs.mLeafBBox)
3367 {
3368 }
3369 
3370 
3371 template<typename SrcTreeT>
3372 void
3374 {
3375  if (threaded) {
3376  tbb::parallel_reduce(mLeafManager.getRange(), *this);
3377  } else {
3378  (*this)(mLeafManager.getRange());
3379  }
3380 }
3381 
3382 
3383 template<typename SrcTreeT>
3384 bool
3385 GenBoundaryMask<SrcTreeT>::neighboringLeaf(const Coord& ijk, const IntTreeAccessorT& acc) const
3386 {
3387  if (acc.probeConstLeaf(ijk)) return true;
3388 
3389  const int dim = IntTreeT::LeafNodeType::DIM;
3390 
3391  // face adjacent neghbours
3392  if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1], ijk[2]))) return true;
3393  if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1], ijk[2]))) return true;
3394  if (acc.probeConstLeaf(Coord(ijk[0], ijk[1] + dim, ijk[2]))) return true;
3395  if (acc.probeConstLeaf(Coord(ijk[0], ijk[1] - dim, ijk[2]))) return true;
3396  if (acc.probeConstLeaf(Coord(ijk[0], ijk[1], ijk[2] + dim))) return true;
3397  if (acc.probeConstLeaf(Coord(ijk[0], ijk[1], ijk[2] - dim))) return true;
3398 
3399  // edge adjacent neighbors
3400  if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1], ijk[2] - dim))) return true;
3401  if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1], ijk[2] - dim))) return true;
3402  if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1], ijk[2] + dim))) return true;
3403  if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1], ijk[2] + dim))) return true;
3404  if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1] + dim, ijk[2]))) return true;
3405  if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1] + dim, ijk[2]))) return true;
3406  if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1] - dim, ijk[2]))) return true;
3407  if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1] - dim, ijk[2]))) return true;
3408  if (acc.probeConstLeaf(Coord(ijk[0], ijk[1] - dim, ijk[2] + dim))) return true;
3409  if (acc.probeConstLeaf(Coord(ijk[0], ijk[1] - dim, ijk[2] - dim))) return true;
3410  if (acc.probeConstLeaf(Coord(ijk[0], ijk[1] + dim, ijk[2] + dim))) return true;
3411  if (acc.probeConstLeaf(Coord(ijk[0], ijk[1] + dim, ijk[2] - dim))) return true;
3412 
3413  // corner adjacent neighbors
3414  if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1] - dim, ijk[2] - dim))) return true;
3415  if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1] - dim, ijk[2] + dim))) return true;
3416  if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1] - dim, ijk[2] + dim))) return true;
3417  if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1] - dim, ijk[2] - dim))) return true;
3418  if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1] + dim, ijk[2] - dim))) return true;
3419  if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1] + dim, ijk[2] + dim))) return true;
3420  if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1] + dim, ijk[2] + dim))) return true;
3421  if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1] + dim, ijk[2] - dim))) return true;
3422 
3423  return false;
3424 }
3425 
3426 
3427 template<typename SrcTreeT>
3428 void
3429 GenBoundaryMask<SrcTreeT>::operator()(const tbb::blocked_range<size_t>& range)
3430 {
3431  Coord ijk;
3432  tree::ValueAccessor<const BoolTreeT> maskAcc(mMaskTree);
3433  tree::ValueAccessor<const IntTreeT> idxAcc(mIdxTree);
3434  tree::ValueAccessor<BoolTreeT> acc(mTree);
3435 
3436  typename SrcTreeT::LeafNodeType::ValueOnCIter iter;
3437 
3438  for (size_t n = range.begin(); n != range.end(); ++n) {
3439 
3440  const typename SrcTreeT::LeafNodeType&
3441  leaf = mLeafManager.leaf(n);
3442 
3443  ijk = leaf.origin();
3444 
3445  if (!mLeafBBox.isInside(ijk) || !neighboringLeaf(ijk, idxAcc)) continue;
3446 
3447  const typename BoolTreeT::LeafNodeType*
3448  maskLeaf = maskAcc.probeConstLeaf(ijk);
3449 
3450  if (!maskLeaf || !leaf.hasSameTopology(maskLeaf)) {
3451  acc.touchLeaf(ijk)->topologyUnion(leaf);
3452  }
3453  }
3454 }
3455 
3456 
3458 
3459 
3460 template<typename TreeT>
3462 {
3463 public:
3464  typedef typename TreeT::template ValueConverter<bool>::Type BoolTreeT;
3465 
3466  typedef typename TreeT::ValueType ValueT;
3467 
3469 
3470  GenTileMask(const std::vector<Vec4i>& tiles, const TreeT& distTree, ValueT iso);
3471 
3472  void run(bool threaded = true);
3473 
3474  BoolTreeT& tree() { return mTree; }
3475 
3477 
3478  GenTileMask(GenTileMask&, tbb::split);
3479  void operator()(const tbb::blocked_range<size_t>&);
3480  void join(GenTileMask& rhs) { mTree.merge(rhs.mTree); }
3481 
3482 private:
3483 
3484  const std::vector<Vec4i>& mTiles;
3485  const TreeT& mDistTree;
3486  ValueT mIsovalue;
3487 
3488  BoolTreeT mTree;
3489 };
3490 
3491 
3492 template<typename TreeT>
3494  const std::vector<Vec4i>& tiles, const TreeT& distTree, ValueT iso)
3495  : mTiles(tiles)
3496  , mDistTree(distTree)
3497  , mIsovalue(iso)
3498  , mTree(false)
3499 {
3500 }
3501 
3502 
3503 template<typename TreeT>
3505  : mTiles(rhs.mTiles)
3506  , mDistTree(rhs.mDistTree)
3507  , mIsovalue(rhs.mIsovalue)
3508  , mTree(false)
3509 {
3510 }
3511 
3512 
3513 template<typename TreeT>
3514 void
3516 {
3517  if (threaded) tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mTiles.size()), *this);
3518  else (*this)(tbb::blocked_range<size_t>(0, mTiles.size()));
3519 }
3520 
3521 
3522 template<typename TreeT>
3523 void
3524 GenTileMask<TreeT>::operator()(const tbb::blocked_range<size_t>& range)
3525 {
3526  tree::ValueAccessor<const TreeT> distAcc(mDistTree);
3527  CoordBBox region, bbox;
3528  Coord ijk, nijk;
3529  bool processRegion = true;
3530  ValueT value;
3531 
3532 
3533  for (size_t n = range.begin(); n != range.end(); ++n) {
3534 
3535  const Vec4i& tile = mTiles[n];
3536 
3537  bbox.min()[0] = tile[0];
3538  bbox.min()[1] = tile[1];
3539  bbox.min()[2] = tile[2];
3540 
3541  bbox.max() = bbox.min();
3542  bbox.max().offset(tile[3]);
3543 
3544  const bool thisInside = (distAcc.getValue(bbox.min()) < mIsovalue);
3545  const int thisDepth = distAcc.getValueDepth(bbox.min());
3546 
3547  // eval x-edges
3548 
3549  ijk = bbox.max();
3550  nijk = ijk;
3551  ++nijk[0];
3552 
3553  processRegion = true;
3554  if (thisDepth >= distAcc.getValueDepth(nijk)) {
3555  processRegion = thisInside != (distAcc.getValue(nijk) < mIsovalue);
3556  }
3557 
3558 
3559  if (processRegion) {
3560  region = bbox;
3561  region.min()[0] = region.max()[0] = ijk[0];
3562  mTree.fill(region, true);
3563  }
3564 
3565 
3566  ijk = bbox.min();
3567  --ijk[0];
3568 
3569  processRegion = true;
3570  if (thisDepth >= distAcc.getValueDepth(ijk)) {
3571  processRegion = !distAcc.probeValue(ijk, value) && thisInside != (value < mIsovalue);
3572  }
3573 
3574  if (processRegion) {
3575  region = bbox;
3576  region.min()[0] = region.max()[0] = ijk[0];
3577  mTree.fill(region, true);
3578  }
3579 
3580 
3581  // eval y-edges
3582 
3583  ijk = bbox.max();
3584  nijk = ijk;
3585  ++nijk[1];
3586 
3587  processRegion = true;
3588  if (thisDepth >= distAcc.getValueDepth(nijk)) {
3589  processRegion = thisInside != (distAcc.getValue(nijk) < mIsovalue);
3590  }
3591 
3592  if (processRegion) {
3593  region = bbox;
3594  region.min()[1] = region.max()[1] = ijk[1];
3595  mTree.fill(region, true);
3596  }
3597 
3598 
3599  ijk = bbox.min();
3600  --ijk[1];
3601 
3602  processRegion = true;
3603  if (thisDepth >= distAcc.getValueDepth(ijk)) {
3604  processRegion = !distAcc.probeValue(ijk, value) && thisInside != (value < mIsovalue);
3605  }
3606 
3607  if (processRegion) {
3608  region = bbox;
3609  region.min()[1] = region.max()[1] = ijk[1];
3610  mTree.fill(region, true);
3611  }
3612 
3613 
3614  // eval z-edges
3615 
3616  ijk = bbox.max();
3617  nijk = ijk;
3618  ++nijk[2];
3619 
3620  processRegion = true;
3621  if (thisDepth >= distAcc.getValueDepth(nijk)) {
3622  processRegion = thisInside != (distAcc.getValue(nijk) < mIsovalue);
3623  }
3624 
3625  if (processRegion) {
3626  region = bbox;
3627  region.min()[2] = region.max()[2] = ijk[2];
3628  mTree.fill(region, true);
3629  }
3630 
3631  ijk = bbox.min();
3632  --ijk[2];
3633 
3634  processRegion = true;
3635  if (thisDepth >= distAcc.getValueDepth(ijk)) {
3636  processRegion = !distAcc.probeValue(ijk, value) && thisInside != (value < mIsovalue);
3637  }
3638 
3639  if (processRegion) {
3640  region = bbox;
3641  region.min()[2] = region.max()[2] = ijk[2];
3642  mTree.fill(region, true);
3643  }
3644 
3645 
3646  ijk = bbox.min();
3647  --ijk[1];
3648  --ijk[2];
3649 
3650  processRegion = true;
3651  if (thisDepth >= distAcc.getValueDepth(ijk)) {
3652  processRegion = !distAcc.probeValue(ijk, value) && thisInside != (value < mIsovalue);
3653  }
3654 
3655  if (processRegion) {
3656  region = bbox;
3657  region.min()[1] = region.max()[1] = ijk[1];
3658  region.min()[2] = region.max()[2] = ijk[2];
3659  mTree.fill(region, true);
3660  }
3661 
3662 
3663  ijk = bbox.min();
3664  --ijk[0];
3665  --ijk[1];
3666 
3667  processRegion = true;
3668  if (thisDepth >= distAcc.getValueDepth(ijk)) {
3669  processRegion = !distAcc.probeValue(ijk, value) && thisInside != (value < mIsovalue);
3670  }
3671 
3672  if (processRegion) {
3673  region = bbox;
3674  region.min()[1] = region.max()[1] = ijk[1];
3675  region.min()[0] = region.max()[0] = ijk[0];
3676  mTree.fill(region, true);
3677  }
3678 
3679  ijk = bbox.min();
3680  --ijk[0];
3681  --ijk[2];
3682 
3683  processRegion = true;
3684  if (thisDepth >= distAcc.getValueDepth(ijk)) {
3685  processRegion = !distAcc.probeValue(ijk, value) && thisInside != (value < mIsovalue);
3686  }
3687 
3688  if (processRegion) {
3689  region = bbox;
3690  region.min()[2] = region.max()[2] = ijk[2];
3691  region.min()[0] = region.max()[0] = ijk[0];
3692  mTree.fill(region, true);
3693  }
3694  }
3695 }
3696 
3697 
3699 
3700 
3701 template<class DistTreeT, class SignTreeT, class IdxTreeT>
3702 inline void
3703 tileData(const DistTreeT& distTree, SignTreeT& signTree, IdxTreeT& idxTree, double iso)
3704 {
3705  typename DistTreeT::ValueOnCIter tileIter(distTree);
3706  tileIter.setMaxDepth(DistTreeT::ValueOnCIter::LEAF_DEPTH - 1);
3707 
3708  if (!tileIter) return; // volume has no active tiles.
3709 
3710  size_t tileCount = 0;
3711  for ( ; tileIter; ++tileIter) {
3712  ++tileCount;
3713  }
3714 
3715  std::vector<Vec4i> tiles(tileCount);
3716 
3717  tileCount = 0;
3718  tileIter = distTree.cbeginValueOn();
3719  tileIter.setMaxDepth(DistTreeT::ValueOnCIter::LEAF_DEPTH - 1);
3720 
3721  CoordBBox bbox;
3722  for (; tileIter; ++tileIter) {
3723  Vec4i& tile = tiles[tileCount++];
3724  tileIter.getBoundingBox(bbox);
3725  tile[0] = bbox.min()[0];
3726  tile[1] = bbox.min()[1];
3727  tile[2] = bbox.min()[2];
3728  tile[3] = bbox.max()[0] - bbox.min()[0];
3729  }
3730 
3731  typename DistTreeT::ValueType isovalue = typename DistTreeT::ValueType(iso);
3732 
3733  GenTileMask<DistTreeT> tileMask(tiles, distTree, isovalue);
3734  tileMask.run();
3735 
3736  typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
3737  typedef tree::LeafManager<BoolTreeT> BoolLeafManagerT;
3738 
3739  BoolLeafManagerT leafs(tileMask.tree());
3740 
3741 
3742  internal::SignData<DistTreeT, BoolLeafManagerT> op(distTree, leafs, isovalue);
3743  op.run();
3744 
3745  signTree.merge(*op.signTree());
3746  idxTree.merge(*op.idxTree());
3747 }
3748 
3749 
3751 
3752 
3753 // Utility class for the volumeToMesh wrapper
3755 {
3756 public:
3757  PointListCopy(const PointList& pointsIn, std::vector<Vec3s>& pointsOut)
3758  : mPointsIn(pointsIn) , mPointsOut(pointsOut)
3759  {
3760  }
3761 
3762  void operator()(const tbb::blocked_range<size_t>& range) const
3763  {
3764  for (size_t n = range.begin(); n < range.end(); ++n) {
3765  mPointsOut[n] = mPointsIn[n];
3766  }
3767  }
3768 
3769 private:
3770  const PointList& mPointsIn;
3771  std::vector<Vec3s>& mPointsOut;
3772 };
3773 
3774 
3775 // Checks if the isovalue is in proximity to the active voxel boundary.
3776 template <typename LeafManagerT>
3777 inline bool
3778 needsActiveVoxePadding(const LeafManagerT& leafs, double iso, double voxelSize)
3779 {
3780  double interiorWidth = 0.0, exteriorWidth = 0.0;
3781  {
3782  typename LeafManagerT::TreeType::LeafNodeType::ValueOffCIter it;
3783  bool foundInterior = false, foundExterior = false;
3784  for (size_t n = 0, N = leafs.leafCount(); n < N; ++n) {
3785 
3786  for (it = leafs.leaf(n).cbeginValueOff(); it; ++it) {
3787  double value = double(it.getValue());
3788  if (value < 0.0) {
3789  interiorWidth = value;
3790  foundInterior = true;
3791  } else if (value > 0.0) {
3792  exteriorWidth = value;
3793  foundExterior = true;
3794  }
3795 
3796  if (foundInterior && foundExterior) break;
3797  }
3798 
3799  if (foundInterior && foundExterior) break;
3800  }
3801 
3802  }
3803 
3804  double minDist = std::min(std::abs(interiorWidth - iso), std::abs(exteriorWidth - iso));
3805  return !(minDist > (2.0 * voxelSize));
3806 }
3807 
3808 
3809 } // end namespace internal
3810 
3811 
3813 
3814 
3815 inline
3817  : mNumQuads(0)
3818  , mNumTriangles(0)
3819  , mQuads(NULL)
3820  , mTriangles(NULL)
3821  , mQuadFlags(NULL)
3822  , mTriangleFlags(NULL)
3823 {
3824 }
3825 
3826 
3827 inline
3828 PolygonPool::PolygonPool(const size_t numQuads, const size_t numTriangles)
3829  : mNumQuads(numQuads)
3830  , mNumTriangles(numTriangles)
3831  , mQuads(new openvdb::Vec4I[mNumQuads])
3832  , mTriangles(new openvdb::Vec3I[mNumTriangles])
3833  , mQuadFlags(new char[mNumQuads])
3834  , mTriangleFlags(new char[mNumTriangles])
3835 {
3836 }
3837 
3838 
3839 inline void
3841 {
3842  resetQuads(rhs.numQuads());
3844 
3845  for (size_t i = 0; i < mNumQuads; ++i) {
3846  mQuads[i] = rhs.mQuads[i];
3847  mQuadFlags[i] = rhs.mQuadFlags[i];
3848  }
3849 
3850  for (size_t i = 0; i < mNumTriangles; ++i) {
3851  mTriangles[i] = rhs.mTriangles[i];
3852  mTriangleFlags[i] = rhs.mTriangleFlags[i];
3853  }
3854 }
3855 
3856 
3857 inline void
3859 {
3860  mNumQuads = size;
3861  mQuads.reset(new openvdb::Vec4I[mNumQuads]);
3862  mQuadFlags.reset(new char[mNumQuads]);
3863 }
3864 
3865 
3866 inline void
3868 {
3869  mNumQuads = 0;
3870  mQuads.reset(NULL);
3871  mQuadFlags.reset(NULL);
3872 }
3873 
3874 
3875 inline void
3877 {
3878  mNumTriangles = size;
3879  mTriangles.reset(new openvdb::Vec3I[mNumTriangles]);
3880  mTriangleFlags.reset(new char[mNumTriangles]);
3881 }
3882 
3883 
3884 inline void
3886 {
3887  mNumTriangles = 0;
3888  mTriangles.reset(NULL);
3889  mTriangleFlags.reset(NULL);
3890 }
3891 
3892 
3893 inline bool
3894 PolygonPool::trimQuads(const size_t n, bool reallocate)
3895 {
3896  if (!(n < mNumQuads)) return false;
3897 
3898  if (reallocate) {
3899 
3900  if (n == 0) {
3901  mQuads.reset(NULL);
3902  } else {
3903 
3904  boost::scoped_array<openvdb::Vec4I> quads(new openvdb::Vec4I[n]);
3905  boost::scoped_array<char> flags(new char[n]);
3906 
3907  for (size_t i = 0; i < n; ++i) {
3908  quads[i] = mQuads[i];
3909  flags[i] = mQuadFlags[i];
3910  }
3911 
3912  mQuads.swap(quads);
3913  mQuadFlags.swap(flags);
3914  }
3915  }
3916 
3917  mNumQuads = n;
3918  return true;
3919 }
3920 
3921 
3922 inline bool
3923 PolygonPool::trimTrinagles(const size_t n, bool reallocate)
3924 {
3925  if (!(n < mNumTriangles)) return false;
3926 
3927  if (reallocate) {
3928 
3929  if (n == 0) {
3930  mTriangles.reset(NULL);
3931  } else {
3932 
3933  boost::scoped_array<openvdb::Vec3I> triangles(new openvdb::Vec3I[n]);
3934  boost::scoped_array<char> flags(new char[n]);
3935 
3936  for (size_t i = 0; i < n; ++i) {
3937  triangles[i] = mTriangles[i];
3938  flags[i] = mTriangleFlags[i];
3939  }
3940 
3941  mTriangles.swap(triangles);
3942  mTriangleFlags.swap(flags);
3943  }
3944  }
3945 
3946  mNumTriangles = n;
3947  return true;
3948 }
3949 
3950 
3952 
3953 
3954 inline VolumeToMesh::VolumeToMesh(double isovalue, double adaptivity)
3955  : mPoints(NULL)
3956  , mPolygons()
3957  , mPointListSize(0)
3958  , mSeamPointListSize(0)
3959  , mPolygonPoolListSize(0)
3960  , mIsovalue(isovalue)
3961  , mPrimAdaptivity(adaptivity)
3962  , mSecAdaptivity(0.0)
3963  , mRefGrid(GridBase::ConstPtr())
3964  , mSurfaceMaskGrid(GridBase::ConstPtr())
3965  , mAdaptivityGrid(GridBase::ConstPtr())
3966  , mAdaptivityMaskTree(TreeBase::ConstPtr())
3967  , mRefSignTree(TreeBase::Ptr())
3968  , mRefIdxTree(TreeBase::Ptr())
3969  , mInvertSurfaceMask(false)
3970  , mPartitions(1)
3971  , mActivePart(0)
3972  , mQuantizedSeamPoints(NULL)
3973  , mPointFlags(0)
3974 {
3975 }
3976 
3977 
3978 inline PointList&
3980 {
3981  return mPoints;
3982 }
3983 
3984 
3985 inline const size_t&
3987 {
3988  return mPointListSize;
3989 }
3990 
3991 
3992 inline PolygonPoolList&
3994 {
3995  return mPolygons;
3996 }
3997 
3998 
3999 inline const PolygonPoolList&
4001 {
4002  return mPolygons;
4003 }
4004 
4005 
4006 inline const size_t&
4008 {
4009  return mPolygonPoolListSize;
4010 }
4011 
4012 
4013 inline void
4014 VolumeToMesh::setRefGrid(const GridBase::ConstPtr& grid, double secAdaptivity)
4015 {
4016  mRefGrid = grid;
4017  mSecAdaptivity = secAdaptivity;
4018 
4019  // Clear out old auxiliary data
4020  mRefSignTree = TreeBase::Ptr();
4021  mRefIdxTree = TreeBase::Ptr();
4022  mSeamPointListSize = 0;
4023  mQuantizedSeamPoints.reset(NULL);
4024 }
4025 
4026 
4027 inline void
4029 {
4030  mSurfaceMaskGrid = mask;
4031  mInvertSurfaceMask = invertMask;
4032 }
4033 
4034 
4035 inline void
4037 {
4038  mAdaptivityGrid = grid;
4039 }
4040 
4041 
4042 inline void
4044 {
4045  mAdaptivityMaskTree = tree;
4046 }
4047 
4048 
4049 inline void
4050 VolumeToMesh::partition(unsigned partitions, unsigned activePart)
4051 {
4052  mPartitions = std::max(partitions, unsigned(1));
4053  mActivePart = std::min(activePart, mPartitions-1);
4054 }
4055 
4056 
4057 inline std::vector<unsigned char>&
4059 {
4060  return mPointFlags;
4061 }
4062 
4063 
4064 inline const std::vector<unsigned char>&
4066 {
4067  return mPointFlags;
4068 }
4069 
4070 
4071 template<typename GridT>
4072 inline void
4073 VolumeToMesh::operator()(const GridT& distGrid)
4074 {
4075  typedef typename GridT::TreeType DistTreeT;
4076  typedef tree::LeafManager<const DistTreeT> DistLeafManagerT;
4077  typedef typename DistTreeT::ValueType DistValueT;
4078 
4079  typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
4080  typedef tree::LeafManager<BoolTreeT> BoolLeafManagerT;
4081  typedef Grid<BoolTreeT> BoolGridT;
4082 
4083  typedef typename DistTreeT::template ValueConverter<Int16>::Type Int16TreeT;
4084  typedef tree::LeafManager<Int16TreeT> Int16LeafManagerT;
4085 
4086  typedef typename DistTreeT::template ValueConverter<int>::Type IntTreeT;
4087  typedef typename DistTreeT::template ValueConverter<float>::Type FloatTreeT;
4088  typedef Grid<FloatTreeT> FloatGridT;
4089 
4090 
4091  const openvdb::math::Transform& transform = distGrid.transform();
4092  const DistTreeT& distTree = distGrid.tree();
4093  const DistValueT isovalue = DistValueT(mIsovalue);
4094 
4095  typename Int16TreeT::Ptr signTreePt;
4096  typename IntTreeT::Ptr idxTreePt;
4097  typename BoolTreeT::Ptr pointMask;
4098 
4099  BoolTreeT valueMask(false), seamMask(false);
4100  const bool adaptive = mPrimAdaptivity > 1e-7 || mSecAdaptivity > 1e-7;
4101  bool maskEdges = false;
4102 
4103 
4104  const BoolGridT * surfaceMask = NULL;
4105  if (mSurfaceMaskGrid && mSurfaceMaskGrid->type() == BoolGridT::gridType()) {
4106  surfaceMask = static_cast<const BoolGridT*>(mSurfaceMaskGrid.get());
4107  }
4108 
4109  const FloatGridT * adaptivityField = NULL;
4110  if (mAdaptivityGrid && mAdaptivityGrid->type() == FloatGridT::gridType()) {
4111  adaptivityField = static_cast<const FloatGridT*>(mAdaptivityGrid.get());
4112  }
4113 
4114  if (mAdaptivityMaskTree && mAdaptivityMaskTree->type() == BoolTreeT::treeType()) {
4115  const BoolTreeT *adaptivityMaskPt =
4116  static_cast<const BoolTreeT*>(mAdaptivityMaskTree.get());
4117  seamMask.topologyUnion(*adaptivityMaskPt);
4118  }
4119 
4120 
4121  // Collect auxiliary data
4122  {
4123  DistLeafManagerT distLeafs(distTree);
4124 
4125  // Check if the isovalue is in proximity to the active voxel boundary.
4126  bool padActiveVoxels = false;
4127  int padVoxels = 3;
4128 
4129  if (distGrid.getGridClass() != GRID_LEVEL_SET) {
4130  padActiveVoxels = true;
4131  } else {
4132  padActiveVoxels = internal::needsActiveVoxePadding(distLeafs,
4133  mIsovalue, transform.voxelSize()[0]);
4134  }
4135 
4136  // always pad the active region for small volumes (the performance hit is neglectable).
4137  if (!padActiveVoxels) {
4138  Coord dim;
4139  distTree.evalActiveVoxelDim(dim);
4140  int maxDim = std::max(std::max(dim[0], dim[1]), dim[2]);
4141  if (maxDim < 1000) {
4142  padActiveVoxels = true;
4143  padVoxels = 1;
4144  }
4145  }
4146 
4147  if (surfaceMask || mPartitions > 1) {
4148 
4149  maskEdges = true;
4150 
4151  if (surfaceMask) {
4152 
4153  { // Mask
4155  *surfaceMask, distLeafs, transform, mInvertSurfaceMask);
4156  masking.run();
4157  valueMask.merge(masking.tree());
4158  }
4159 
4160  if (mPartitions > 1) { // Partition
4161  tree::LeafManager<BoolTreeT> leafs(valueMask);
4162  leafs.foreach(internal::PartOp(leafs.leafCount() , mPartitions, mActivePart));
4163  valueMask.pruneInactive();
4164  }
4165 
4166  } else { // Partition
4167 
4168  internal::PartGen<DistTreeT> partitioner(distLeafs, mPartitions, mActivePart);
4169  partitioner.run();
4170  valueMask.merge(partitioner.tree());
4171  }
4172 
4173  {
4174  if (padActiveVoxels) tools::dilateVoxels(valueMask, padVoxels);
4175  BoolLeafManagerT leafs(valueMask);
4176 
4178  signDataOp(distTree, leafs, isovalue);
4179  signDataOp.run();
4180 
4181  signTreePt = signDataOp.signTree();
4182  idxTreePt = signDataOp.idxTree();
4183  }
4184 
4185  {
4186  internal::GenBoundaryMask<DistTreeT> boundary(distLeafs, valueMask, *idxTreePt);
4187  boundary.run();
4188 
4189  BoolLeafManagerT bleafs(boundary.tree());
4190 
4192  signDataOp(distTree, bleafs, isovalue);
4193  signDataOp.run();
4194 
4195  signTreePt->merge(*signDataOp.signTree());
4196  idxTreePt->merge(*signDataOp.idxTree());
4197  }
4198 
4199  } else {
4200 
4201  // Collect voxel-sign configurations
4202  if (padActiveVoxels) {
4203 
4204  BoolTreeT regionMask(false);
4205  regionMask.topologyUnion(distTree);
4206  tools::dilateVoxels(regionMask, padVoxels);
4207 
4208  BoolLeafManagerT leafs(regionMask);
4209 
4211  signDataOp(distTree, leafs, isovalue);
4212  signDataOp.run();
4213 
4214  signTreePt = signDataOp.signTree();
4215  idxTreePt = signDataOp.idxTree();
4216  } else {
4217 
4219  signDataOp(distTree, distLeafs, isovalue);
4220  signDataOp.run();
4221 
4222  signTreePt = signDataOp.signTree();
4223  idxTreePt = signDataOp.idxTree();
4224  }
4225  }
4226 
4227  }
4228 
4229 
4230  // Collect auxiliary data from active tiles
4231  internal::tileData(distTree, *signTreePt, *idxTreePt, isovalue);
4232 
4233  // Optionally collect auxiliary data from a reference level set.
4234  Int16TreeT *refSignTreePt = NULL;
4235  IntTreeT *refIdxTreePt = NULL;
4236  const DistTreeT *refDistTreePt = NULL;
4237 
4238  if (mRefGrid && mRefGrid->type() == GridT::gridType()) {
4239 
4240  const GridT* refGrid = static_cast<const GridT*>(mRefGrid.get());
4241  refDistTreePt = &refGrid->tree();
4242 
4243  // Collect and cache auxiliary data from the reference grid.
4244  if (!mRefSignTree && !mRefIdxTree) {
4245 
4246  DistLeafManagerT refDistLeafs(*refDistTreePt);
4248  signDataOp(*refDistTreePt, refDistLeafs, isovalue);
4249 
4250  signDataOp.run();
4251 
4252  mRefSignTree = signDataOp.signTree();
4253  mRefIdxTree = signDataOp.idxTree();
4254  }
4255 
4256  // Get cached auxiliary data
4257  if (mRefSignTree && mRefIdxTree) {
4258  refSignTreePt = static_cast<Int16TreeT*>(mRefSignTree.get());
4259  refIdxTreePt = static_cast<IntTreeT*>(mRefIdxTree.get());
4260  }
4261  }
4262 
4263 
4264  // Process auxiliary data
4265  Int16LeafManagerT signLeafs(*signTreePt);
4266 
4267  if (maskEdges) {
4268  signLeafs.foreach(internal::MaskEdges<BoolTreeT>(valueMask));
4269  valueMask.clear();
4270  }
4271 
4272 
4273  // Generate the seamline mask
4274  if (refSignTreePt) {
4275  internal::GenSeamMask<Int16TreeT, Int16LeafManagerT> seamOp(signLeafs, *refSignTreePt);
4276  seamOp.run();
4277 
4278  tools::dilateVoxels(seamOp.mask(), 3);
4279  signLeafs.foreach(internal::TagSeamEdges<BoolTreeT>(seamOp.mask()));
4280 
4281  seamMask.merge(seamOp.mask());
4282  }
4283 
4284 
4285  std::vector<size_t> regions(signLeafs.leafCount(), 0);
4286  if (regions.empty()) return;
4287 
4288  if (adaptive) {
4289 
4291  signLeafs, *signTreePt, distTree, *idxTreePt, isovalue, DistValueT(mPrimAdaptivity));
4292 
4293  if (adaptivityField) {
4294  merge.setSpatialAdaptivity(transform, *adaptivityField);
4295  }
4296 
4297  if (refSignTreePt || mAdaptivityMaskTree) {
4298  merge.setAdaptivityMask(&seamMask);
4299  }
4300 
4301  if (refSignTreePt) {
4302  merge.setRefData(refSignTreePt, DistValueT(mSecAdaptivity));
4303  }
4304 
4305  merge.run();
4306 
4307  signLeafs.foreach(internal::CountRegions<IntTreeT>(*idxTreePt, regions));
4308 
4309  } else {
4310  signLeafs.foreach(internal::CountPoints(regions));
4311  }
4312 
4313 
4314  {
4315  mPointListSize = 0;
4316  size_t tmp = 0;
4317  for (size_t n = 0, N = regions.size(); n < N; ++n) {
4318  tmp = regions[n];
4319  regions[n] = mPointListSize;
4320  mPointListSize += tmp;
4321  }
4322  }
4323 
4324 
4325  // Generate the unique point list
4326  mPoints.reset(new openvdb::Vec3s[mPointListSize]);
4327  mPointFlags.clear();
4328 
4329  // Generate seam line sample points
4330  if (refSignTreePt && refIdxTreePt) {
4331 
4332  if (mSeamPointListSize == 0) {
4333 
4334  std::vector<size_t> pointMap;
4335 
4336  {
4337  Int16LeafManagerT refSignLeafs(*refSignTreePt);
4338  pointMap.resize(refSignLeafs.leafCount(), 0);
4339 
4340  refSignLeafs.foreach(internal::CountPoints(pointMap));
4341 
4342  size_t tmp = 0;
4343  for (size_t n = 0, N = pointMap.size(); n < N; ++n) {
4344  tmp = pointMap[n];
4345  pointMap[n] = mSeamPointListSize;
4346  mSeamPointListSize += tmp;
4347  }
4348  }
4349 
4350  if (!pointMap.empty() && mSeamPointListSize != 0) {
4351 
4352  mQuantizedSeamPoints.reset(new uint32_t[mSeamPointListSize]);
4353  memset(mQuantizedSeamPoints.get(), 0, sizeof(uint32_t) * mSeamPointListSize);
4354 
4355  typedef tree::LeafManager<IntTreeT> IntLeafManagerT;
4356 
4357  IntLeafManagerT refIdxLeafs(*refIdxTreePt);
4358  refIdxLeafs.foreach(internal::MapPoints<Int16TreeT>(pointMap, *refSignTreePt));
4359  }
4360  }
4361 
4362  if (mSeamPointListSize != 0) {
4363  signLeafs.foreach(internal::SeamWeights<DistTreeT>(
4364  distTree, *refSignTreePt, *refIdxTreePt, mQuantizedSeamPoints, mIsovalue));
4365  }
4366  }
4367 
4368 
4370  pointOp(signLeafs, distTree, *idxTreePt, mPoints, regions, transform, mIsovalue);
4371 
4372 
4373  if (mSeamPointListSize != 0) {
4374  mPointFlags.resize(mPointListSize);
4375  pointOp.setRefData(refSignTreePt, refDistTreePt, refIdxTreePt,
4376  &mQuantizedSeamPoints, &mPointFlags);
4377  }
4378 
4379  pointOp.run();
4380 
4381 
4382  mPolygonPoolListSize = signLeafs.leafCount();
4383  mPolygons.reset(new PolygonPool[mPolygonPoolListSize]);
4384 
4385 
4386  if (adaptive) {
4387 
4389  mesher(signLeafs, *signTreePt, *idxTreePt, mPolygons, Index32(mPointListSize));
4390 
4391  mesher.setRefSignTree(refSignTreePt);
4392  mesher.run();
4393 
4394  } else {
4395 
4397  mesher(signLeafs, *signTreePt, *idxTreePt, mPolygons, Index32(mPointListSize));
4398 
4399  mesher.setRefSignTree(refSignTreePt);
4400  mesher.run();
4401  }
4402 
4403  // Clean up unused points, only necessary if masking and/or
4404  // automatic mesh partitioning is enabled.
4405  if ((surfaceMask || mPartitions > 1) && mPointListSize > 0) {
4406 
4407  // Flag used points
4408  std::vector<unsigned char> usedPointMask(mPointListSize, 0);
4409 
4410  internal::FlagUsedPoints flagPoints(mPolygons, mPolygonPoolListSize, usedPointMask);
4411  flagPoints.run();
4412 
4413  // Create index map
4414  std::vector<unsigned> indexMap(mPointListSize);
4415  size_t usedPointCount = 0;
4416  for (size_t p = 0; p < mPointListSize; ++p) {
4417  if (usedPointMask[p]) indexMap[p] = usedPointCount++;
4418  }
4419 
4420  if (usedPointCount < mPointListSize) {
4421 
4422  // move points
4424  newPointList(new openvdb::Vec3s[usedPointCount]);
4425 
4426  internal::MovePoints movePoints(newPointList, mPoints, indexMap, usedPointMask);
4427  movePoints.run();
4428 
4429  mPointListSize = usedPointCount;
4430  mPoints.reset(newPointList.release());
4431 
4432  // update primitives
4433  internal::RemapIndices remap(mPolygons, mPolygonPoolListSize, indexMap);
4434  remap.run();
4435  }
4436  }
4437 
4438 
4439  // Subdivide nonplanar quads near the seamline edges
4440  // todo: thread and clean up
4441  if (refSignTreePt || refIdxTreePt || refDistTreePt) {
4442  std::vector<Vec3s> newPoints;
4443 
4444  for (size_t n = 0; n < mPolygonPoolListSize; ++n) {
4445 
4446  PolygonPool& polygons = mPolygons[n];
4447 
4448  std::vector<size_t> nonPlanarQuads;
4449  nonPlanarQuads.reserve(polygons.numQuads());
4450 
4451  for (size_t i = 0; i < polygons.numQuads(); ++i) {
4452 
4453  char& flags = polygons.quadFlags(i);
4454 
4455  if ((flags & POLYFLAG_FRACTURE_SEAM) && !(flags & POLYFLAG_EXTERIOR)) {
4456 
4457  openvdb::Vec4I& quad = polygons.quad(i);
4458 
4459  const bool edgePoly = mPointFlags[quad[0]] || mPointFlags[quad[1]]
4460  || mPointFlags[quad[2]] || mPointFlags[quad[3]];
4461 
4462  if (!edgePoly) continue;
4463 
4464  const Vec3s& p0 = mPoints[quad[0]];
4465  const Vec3s& p1 = mPoints[quad[1]];
4466  const Vec3s& p2 = mPoints[quad[2]];
4467  const Vec3s& p3 = mPoints[quad[3]];
4468 
4469  if (!internal::isPlanarQuad(p0, p1, p2, p3, 1e-6f)) {
4470  nonPlanarQuads.push_back(i);
4471  }
4472  }
4473  }
4474 
4475 
4476  if (!nonPlanarQuads.empty()) {
4477 
4478  PolygonPool tmpPolygons;
4479 
4480  tmpPolygons.resetQuads(polygons.numQuads() - nonPlanarQuads.size());
4481  tmpPolygons.resetTriangles(polygons.numTriangles() + 4 * nonPlanarQuads.size());
4482 
4483  size_t triangleIdx = 0;
4484  for (size_t i = 0; i < nonPlanarQuads.size(); ++i) {
4485 
4486  size_t& quadIdx = nonPlanarQuads[i];
4487 
4488  openvdb::Vec4I& quad = polygons.quad(quadIdx);
4489  char& quadFlags = polygons.quadFlags(quadIdx);
4490  //quadFlags |= POLYFLAG_SUBDIVIDED;
4491 
4492  Vec3s centroid = (mPoints[quad[0]] + mPoints[quad[1]] +
4493  mPoints[quad[2]] + mPoints[quad[3]]) * 0.25;
4494 
4495  size_t pointIdx = newPoints.size() + mPointListSize;
4496 
4497  newPoints.push_back(centroid);
4498 
4499 
4500  {
4501  Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
4502 
4503  triangle[0] = quad[0];
4504  triangle[1] = pointIdx;
4505  triangle[2] = quad[3];
4506 
4507  tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
4508 
4509  if (mPointFlags[triangle[0]] || mPointFlags[triangle[2]]) {
4510  tmpPolygons.triangleFlags(triangleIdx) |= POLYFLAG_SUBDIVIDED;
4511  }
4512  }
4513 
4514  ++triangleIdx;
4515 
4516  {
4517  Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
4518 
4519  triangle[0] = quad[0];
4520  triangle[1] = quad[1];
4521  triangle[2] = pointIdx;
4522 
4523  tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
4524 
4525  if (mPointFlags[triangle[0]] || mPointFlags[triangle[1]]) {
4526  tmpPolygons.triangleFlags(triangleIdx) |= POLYFLAG_SUBDIVIDED;
4527  }
4528  }
4529 
4530  ++triangleIdx;
4531 
4532  {
4533  Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
4534 
4535  triangle[0] = quad[1];
4536  triangle[1] = quad[2];
4537  triangle[2] = pointIdx;
4538 
4539  tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
4540 
4541  if (mPointFlags[triangle[0]] || mPointFlags[triangle[1]]) {
4542  tmpPolygons.triangleFlags(triangleIdx) |= POLYFLAG_SUBDIVIDED;
4543  }
4544  }
4545 
4546 
4547  ++triangleIdx;
4548 
4549  {
4550  Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
4551 
4552  triangle[0] = quad[2];
4553  triangle[1] = quad[3];
4554  triangle[2] = pointIdx;
4555 
4556  tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
4557 
4558  if (mPointFlags[triangle[0]] || mPointFlags[triangle[1]]) {
4559  tmpPolygons.triangleFlags(triangleIdx) |= POLYFLAG_SUBDIVIDED;
4560  }
4561  }
4562 
4563  ++triangleIdx;
4564 
4565  quad[0] = util::INVALID_IDX;
4566  }
4567 
4568 
4569  for (size_t i = 0; i < polygons.numTriangles(); ++i) {
4570  tmpPolygons.triangle(triangleIdx) = polygons.triangle(i);
4571  tmpPolygons.triangleFlags(triangleIdx) = polygons.triangleFlags(i);
4572  ++triangleIdx;
4573  }
4574 
4575 
4576  size_t quadIdx = 0;
4577  for (size_t i = 0; i < polygons.numQuads(); ++i) {
4578  openvdb::Vec4I& quad = polygons.quad(i);
4579 
4580  if (quad[0] != util::INVALID_IDX) {
4581  tmpPolygons.quad(quadIdx) = quad;
4582  tmpPolygons.quadFlags(quadIdx) = polygons.quadFlags(i);
4583  ++quadIdx;
4584  }
4585  }
4586 
4587 
4588  polygons.copy(tmpPolygons);
4589  }
4590 
4591  }
4592 
4593 
4594  if (!newPoints.empty()) {
4595 
4596  size_t newPointCount = newPoints.size() + mPointListSize;
4597 
4599  newPointList(new openvdb::Vec3s[newPointCount]);
4600 
4601  for (size_t i = 0; i < mPointListSize; ++i) {
4602  newPointList.get()[i] = mPoints[i];
4603  }
4604 
4605  for (size_t i = mPointListSize; i < newPointCount; ++i) {
4606  newPointList.get()[i] = newPoints[i - mPointListSize];
4607  }
4608 
4609  mPointListSize = newPointCount;
4610  mPoints.reset(newPointList.release());
4611  mPointFlags.resize(mPointListSize, 0);
4612  }
4613  }
4614 }
4615 
4616 
4618 
4619 
4620 template<typename GridType>
4621 inline void
4623  const GridType& grid,
4624  std::vector<Vec3s>& points,
4625  std::vector<Vec3I>& triangles,
4626  std::vector<Vec4I>& quads,
4627  double isovalue,
4628  double adaptivity)
4629 {
4630  VolumeToMesh mesher(isovalue, adaptivity);
4631  mesher(grid);
4632 
4633  // Preallocate the point list
4634  points.clear();
4635  points.resize(mesher.pointListSize());
4636 
4637  { // Copy points
4638  internal::PointListCopy ptnCpy(mesher.pointList(), points);
4639  tbb::parallel_for(tbb::blocked_range<size_t>(0, points.size()), ptnCpy);
4640  mesher.pointList().reset(NULL);
4641  }
4642 
4643  PolygonPoolList& polygonPoolList = mesher.polygonPoolList();
4644 
4645  { // Preallocate primitive lists
4646  size_t numQuads = 0, numTriangles = 0;
4647  for (size_t n = 0, N = mesher.polygonPoolListSize(); n < N; ++n) {
4648  openvdb::tools::PolygonPool& polygons = polygonPoolList[n];
4649  numTriangles += polygons.numTriangles();
4650  numQuads += polygons.numQuads();
4651  }
4652 
4653  triangles.clear();
4654  triangles.resize(numTriangles);
4655  quads.clear();
4656  quads.resize(numQuads);
4657  }
4658 
4659  // Copy primitives
4660  size_t qIdx = 0, tIdx = 0;
4661  for (size_t n = 0, N = mesher.polygonPoolListSize(); n < N; ++n) {
4662  openvdb::tools::PolygonPool& polygons = polygonPoolList[n];
4663 
4664  for (size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4665  quads[qIdx++] = polygons.quad(i);
4666  }
4667 
4668  for (size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
4669  triangles[tIdx++] = polygons.triangle(i);
4670  }
4671  }
4672 }
4673 
4674 
4675 template<typename GridType>
4676 void
4678  const GridType& grid,
4679  std::vector<Vec3s>& points,
4680  std::vector<Vec4I>& quads,
4681  double isovalue)
4682 {
4683  std::vector<Vec3I> triangles(0);
4684  volumeToMesh(grid,points, triangles, quads, isovalue, 0.0);
4685 }
4686 
4687 
4689 
4690 
4691 } // namespace tools
4692 } // namespace OPENVDB_VERSION_NAME
4693 } // namespace openvdb
4694 
4695 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
4696 
4697 // Copyright (c) 2012-2013 DreamWorks Animation LLC
4698 // All rights reserved. This software is distributed under the
4699 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
void setRefSignTree(const Int16TreeT *r)
Definition: VolumeToMesh.h:2633
OPENVDB_API Hermite min(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
Definition: VolumeToMesh.h:2761
bool trimQuads(const size_t n, bool reallocate=false)
Definition: VolumeToMesh.h:3894
bool trimTrinagles(const size_t n, bool reallocate=false)
Definition: VolumeToMesh.h:3923
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:4014
const size_t & numTriangles() const
Definition: VolumeToMesh.h:133
void operator()(const tbb::blocked_range< size_t > &)
Definition: VolumeToMesh.h:3524
tree::ValueAccessor< const TreeT > AccessorT
Definition: VolumeToMesh.h:1907
const char & quadFlags(size_t n) const
Definition: VolumeToMesh.h:142
int16_t Int16
Definition: Types.h:58
Definition: VolumeToMesh.h:894
void correctCellSigns(unsigned char &signs, unsigned char 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:657
VolumeToMesh(double isovalue=0, double adaptivity=0)
Definition: VolumeToMesh.h:3954
PolygonPool()
Definition: VolumeToMesh.h:3816
tree::ValueAccessor< const Int16TreeT > Int16AccessorT
Definition: VolumeToMesh.h:2029
CountRegions(IntTreeT &idxTree, std::vector< size_t > &regions)
Definition: VolumeToMesh.h:1117
MovePoints(internal::UniquePtr< openvdb::Vec3s >::type &newPointList, const PointList &oldPointList, const std::vector< unsigned > &indexMap, const std::vector< unsigned char > &usedPointMask)
Definition: VolumeToMesh.h:3135
boost::scoped_array< uint32_t > QuantizedPointList
Definition: VolumeToMesh.h:1915
void setAdaptivityMask(const BoolTreeT *mask)
Definition: VolumeToMesh.h:2116
TagSeamEdges(const TreeT &tree)
Definition: VolumeToMesh.h:2964
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:328
MaskEdges(const BoolTreeT &valueMask)
Definition: VolumeToMesh.h:2995
void operator()(const tbb::blocked_range< size_t > &) const
Definition: VolumeToMesh.h:2674
boost::scoped_array< openvdb::Vec3s > PointList
Point and primitive list types.
Definition: VolumeToMesh.h:167
TreeT::ValueType ValueT
Definition: VolumeToMesh.h:897
LeafManagerT::TreeType::template ValueConverter< Int16 >::Type Int16TreeT
Definition: VolumeToMesh.h:2619
void operator()(const tbb::blocked_range< size_t > &)
Definition: VolumeToMesh.h:982
BoolTreeT & tree()
Definition: VolumeToMesh.h:3198
bool isMergable(LeafType &leaf, const Coord &start, int dim, typename LeafType::ValueType::value_type adaptivity)
Definition: VolumeToMesh.h:856
double evalRoot(double v0, double v1, double iso)
Definition: VolumeToMesh.h:1165
const unsigned char sAmbiguousFace[256]
Contains the ambiguous face index for certain cell configuration.
Definition: VolumeToMesh.h:413
bool isPlanarQuad(const Vec3d &p0, const Vec3d &p1, const Vec3d &p2, const Vec3d &p3, double epsilon=0.001)
Definition: VolumeToMesh.h:519
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: VolumeToMesh.h:3762
boost::shared_ptr< const TreeBase > ConstPtr
Definition: Tree.h:68
Grid< BoolTreeT > BoolGridT
Definition: VolumeToMesh.h:3187
void addPrim(const Vec4I &verts, bool reverse, char flags=0)
Definition: VolumeToMesh.h:2382
void mergeVoxels(LeafType &leaf, const Coord &start, int dim, int regionId)
Definition: VolumeToMesh.h:835
std::auto_ptr< T > type
Definition: VolumeToMesh.h:390
void setSurfaceMask(const GridBase::ConstPtr &mask, bool invertMask=false)
Definition: VolumeToMesh.h:4028
void partition(unsigned partitions=1, unsigned activePart=0)
Subdivide volume and mesh into disjoint parts.
Definition: VolumeToMesh.h:4050
openvdb::Vec3I & triangle(size_t n)
Definition: VolumeToMesh.h:135
RemapIndices(PolygonPoolList &polygons, size_t polyListCount, const std::vector< unsigned > &indexMap)
Definition: VolumeToMesh.h:3082
MapPoints(std::vector< size_t > &pointList, const Int16TreeT &signTree)
Definition: VolumeToMesh.h:1081
TreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: VolumeToMesh.h:2872
math::Transform & transform()
Return a reference to this grid's transform, which might be shared with other grids.
Definition: Grid.h:321
SignData(const TreeT &distTree, const LeafManagerT &leafs, ValueT iso)
Definition: VolumeToMesh.h:944
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:4677
TreeT::template ValueConverter< int >::Type IntTreeT
Definition: VolumeToMesh.h:2023
Computes the point list indices for the index tree.
Definition: VolumeToMesh.h:1076
Definition: VolumeToMesh.h:396
T & z()
Definition: Vec3.h:96
void join(const SignData &rhs)
Definition: VolumeToMesh.h:920
void init(const size_t upperBound, PolygonPool &quadPool)
Definition: VolumeToMesh.h:2334
LeafManagerT::TreeType::template ValueConverter< int >::Type IntTreeT
Definition: VolumeToMesh.h:2618
Index32 Index
Definition: Types.h:57
void join(GenSeamMask &rhs)
Definition: VolumeToMesh.h:2886
LeafNodeT * probeLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains voxel (x, y, z), or NULL if no such node exists...
Definition: ValueAccessor.h:378
void operator()(const tbb::blocked_range< size_t > &)
Definition: VolumeToMesh.h:2926
bool normalize(T eps=T(1.0e-7))
this = normalized this
Definition: Vec3.h:344
Definition: VolumeToMesh.h:396
OPENVDB_API const Index32 INVALID_IDX
boost::scoped_array< PolygonPool > PolygonPoolList
Point and primitive list types.
Definition: VolumeToMesh.h:168
Definition: Mat.h:146
Vec3d unpackPoint(uint32_t data)
Utility methods for point quantization.
Definition: VolumeToMesh.h:573
Vec3d worldToIndex(const Vec3d &xyz) const
Apply this transformation to the given coordinates.
Definition: Transform.h:136
Mat3< double > Mat3d
Definition: Mat3.h:666
const openvdb::Vec3I & triangle(size_t n) const
Definition: VolumeToMesh.h:136
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:1413
Definition: VolumeToMesh.h:397
Vec4< int32_t > Vec4i
Definition: Vec4.h:543
unsigned char evalCellSigns(const AccessorT &accessor, const Coord &ijk, typename AccessorT::ValueType iso)
General method that computes the cell-sign configuration at the given ijk coordinate.
Definition: VolumeToMesh.h:594
void resetTriangles(size_t size)
Definition: VolumeToMesh.h:3876
GenTileMask(const std::vector< Vec4i > &tiles, const TreeT &distTree, ValueT iso)
Definition: VolumeToMesh.h:3493
bool needsActiveVoxePadding(const LeafManagerT &leafs, double iso, double voxelSize)
Definition: VolumeToMesh.h:3778
void setRefData(const Int16TreeT *signTree, ValueT adaptivity)
Definition: VolumeToMesh.h:2123
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:1309
Mat3 transpose() const
returns transpose of this
Definition: Mat3.h:455
void operator()(LeafNodeType &leaf, size_t) const
Definition: VolumeToMesh.h:2967
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: VolumeToMesh.h:3047
void tileData(const DistTreeT &distTree, SignTreeT &signTree, IdxTreeT &idxTree, double iso)
Definition: VolumeToMesh.h:3703
Grid< FloatTreeT > FloatGridT
Definition: VolumeToMesh.h:2032
void run(bool threaded=true)
Definition: VolumeToMesh.h:974
tree::ValueAccessor< Int16TreeT > Int16AccessorT
Definition: VolumeToMesh.h:904
SrcTreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: VolumeToMesh.h:3184
const bool sAdaptable[256]
Used to quickly determine if a given cell is adaptable.
Definition: VolumeToMesh.h:401
IntTreeT::Ptr idxTree() const
Definition: VolumeToMesh.h:914
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: VolumeToMesh.h:3158
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:1584
tree::ValueAccessor< const TreeT > AccessorT
Definition: VolumeToMesh.h:2021
tree::LeafManager< const SrcTreeT > LeafManagerT
Definition: VolumeToMesh.h:3183
TreeT::template ValueConverter< Int16 >::Type Int16TreeT
Definition: VolumeToMesh.h:1912
GenPolygons(const LeafManagerT &signLeafs, const Int16TreeT &signTree, const IntTreeT &idxTree, PolygonPoolList &polygons, Index32 pointListSize)
Definition: VolumeToMesh.h:2652
void operator()(const tbb::blocked_range< size_t > &) const
Definition: VolumeToMesh.h:2132
IntTreeT::LeafNodeType IntLeafT
Definition: VolumeToMesh.h:1115
void clearQuads()
Definition: VolumeToMesh.h:3867
tree::ValueAccessor< const TreeT > AccessorT
Definition: VolumeToMesh.h:898
tree::ValueAccessor< const TreeT > AccessorT
Definition: VolumeToMesh.h:1637
Definition: VolumeToMesh.h:2869
void done()
Definition: VolumeToMesh.h:2420
tree::ValueAccessor< const Int16TreeT > Int16AccessorT
Definition: VolumeToMesh.h:2622
char & triangleFlags(size_t n)
Definition: VolumeToMesh.h:144
Counts the total number of points per leaf, accounts for cells with multiple points.
Definition: VolumeToMesh.h:1051
tree::ValueAccessor< IntTreeT > IntAccessorT
Definition: VolumeToMesh.h:2024
Abstract base class for typed grids.
Definition: Grid.h:103
GenBoundaryMask(const LeafManagerT &leafs, const BoolTreeT &, const IntTreeT &)
Definition: VolumeToMesh.h:3348
TreeT::template ValueConverter< int >::Type IntTreeT
Definition: VolumeToMesh.h:1909
TreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: VolumeToMesh.h:3464
Definition: Types.h:170
tree::ValueAccessor< const Int16TreeT > Int16CAccessorT
Definition: VolumeToMesh.h:1644
Definition: VolumeToMesh.h:104
void join(GenBoundaryMask &rhs)
Definition: VolumeToMesh.h:3331
void join(GenTopologyMask &rhs)
Definition: VolumeToMesh.h:3207
tree::ValueAccessor< BoolTreeT > BoolAccessorT
Definition: VolumeToMesh.h:3186
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: ValueAccessor.h:210
Definition: VolumeToMesh.h:397
Definition: VolumeToMesh.h:2786
bool isNonManifold(const AccessorT &accessor, const Coord &ijk, typename AccessorT::ValueType isovalue, const int dim)
Definition: VolumeToMesh.h:684
TreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: VolumeToMesh.h:2026
void run(bool threaded=true)
Definition: VolumeToMesh.h:3515
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:1217
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:54
const openvdb::Vec4I & quad(size_t n) const
Definition: VolumeToMesh.h:130
void join(GenTileMask &rhs)
Definition: VolumeToMesh.h:3480
#define OPENVDB_VERSION_NAME
Definition: version.h:45
Int16TreeT::Ptr signTree() const
Definition: VolumeToMesh.h:913
PartOp(size_t leafCount, size_t partitions, size_t activePart)
Definition: VolumeToMesh.h:2764
const size_t & pointListSize() const
Definition: VolumeToMesh.h:3986
T dot(const Vec3< T > &v) const
Dot product.
Definition: Vec3.h:199
void operator()(const GridT &)
Main call.
Definition: VolumeToMesh.h:4073
void addPrim(const Vec4I &verts, bool reverse, char flags=0)
Definition: VolumeToMesh.h:2341
bool probeValue(const Coord &xyz, ValueType &value) const
Return the active state of the voxel as well as its value.
Definition: ValueAccessor.h:220
boost::shared_ptr< const GridBase > ConstPtr
Definition: Grid.h:107
Gradient operators defined in index space of various orders.
Definition: Operators.h:122
Definition: VolumeToMesh.h:385
char & quadFlags(size_t n)
Definition: VolumeToMesh.h:141
void operator()(const tbb::blocked_range< size_t > &)
Definition: VolumeToMesh.h:3429
const size_t & polygonPoolListSize() const
Definition: VolumeToMesh.h:4007
openvdb::Vec4I & quad(size_t n)
Definition: VolumeToMesh.h:129
OPENVDB_API Hermite max(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
Collection of quads and triangles.
Definition: VolumeToMesh.h:108
Coord nearestCoord(const Vec3d &voxelCoord)
Return voxelCoord rounded to the closest integer coordinates.
Definition: util/Util.h:55
void operator()(LeafNodeType &leaf, size_t leafIndex) const
Definition: VolumeToMesh.h:2772
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:383
tree::ValueAccessor< const Int16TreeT > Int16AccessorT
Definition: VolumeToMesh.h:1079
uint32_t packPoint(const Vec3d &v)
Utility methods for point quantization.
Definition: VolumeToMesh.h:557
size_t leafCount() const
Return the number of leaf nodes.
Definition: LeafManager.h:295
Base class for typed trees.
Definition: Tree.h:64
BoolTreeT & tree()
Definition: VolumeToMesh.h:3325
void init(const size_t upperBound, PolygonPool &polygonPool)
Definition: VolumeToMesh.h:2372
void copy(const PolygonPool &rhs)
Definition: VolumeToMesh.h:3840
void run(bool threaded=true)
Definition: VolumeToMesh.h:3147
Mesh any scalar grid that has a continuous isosurface.
Definition: VolumeToMesh.h:176
void operator()(const tbb::blocked_range< size_t > &)
Definition: VolumeToMesh.h:2849
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:427
void operator()(LeafNodeType &leaf, size_t leafIndex) const
Definition: VolumeToMesh.h:1057
FlagUsedPoints(const PolygonPoolList &polygons, size_t polyListCount, std::vector< unsigned char > &usedPointMask)
Definition: VolumeToMesh.h:3028
GenTopologyMask(const BoolGridT &mask, const LeafManagerT &srcLeafs, const math::Transform &srcXForm, bool invertMask)
Definition: VolumeToMesh.h:3220
LeafManagerT::TreeType::template ValueConverter< Int16 >::Type Int16TreeT
Definition: VolumeToMesh.h:2028
std::vector< unsigned char > & pointFlags()
Definition: VolumeToMesh.h:4058
Definition: VolumeToMesh.h:2615
const size_t & numQuads() const
Definition: VolumeToMesh.h:127
void run(bool threaded=true)
Definition: VolumeToMesh.h:3244
tree::LeafManager< const SrcTreeT > LeafManagerT
Definition: VolumeToMesh.h:3317
boost::shared_ptr< TreeBase > Ptr
Definition: Tree.h:67
tree::ValueAccessor< IntTreeT > IntAccessorT
Definition: VolumeToMesh.h:1640
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:1571
uint32_t Index32
Definition: Types.h:55
tree::ValueAccessor< BoolTreeT > BoolAccessorT
Definition: VolumeToMesh.h:2791
void collectCornerValues(const LeafT &leaf, const Index offset, std::vector< double > &values)
Extracts the eight corner values for leaf inclusive cells.
Definition: VolumeToMesh.h:1171
void resetQuads(size_t size)
Definition: VolumeToMesh.h:3858
void operator()(const tbb::blocked_range< size_t > &)
Definition: VolumeToMesh.h:3256
void foreach(const LeafOp &op, bool threaded=true, size_t grainSize=1)
Threaded method that applies a user-supplied functor to each leaf node in the LeafManager.
Definition: LeafManager.h:463
TreeT::ValueType ValueT
Definition: VolumeToMesh.h:3466
Vec3< double > Vec3d
Definition: Vec3.h:625
unsigned char 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:621
void run(bool threaded=true)
Definition: VolumeToMesh.h:3036
Vec3< float > Vec3s
Definition: Vec3.h:624
tree::ValueAccessor< IntTreeT > IntAccessorT
Definition: VolumeToMesh.h:1114
void done()
Definition: VolumeToMesh.h:2356
void run(bool threaded=true)
Definition: VolumeToMesh.h:3373
SrcTreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: VolumeToMesh.h:2790
void setAdaptivityMask(const TreeBase::ConstPtr &tree)
Definition: VolumeToMesh.h:4043
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)
Adaptively mesh any scalar grid that has a continuous isosurface.
Definition: VolumeToMesh.h:4622
void run(bool threaded=true)
Definition: VolumeToMesh.h:2097
GenSeamMask(const LeafManagerT &leafs, const TreeT &tree)
Definition: VolumeToMesh.h:2898
TreeT::template ValueConverter< int >::Type IntTreeT
Definition: VolumeToMesh.h:1639
void operator()(LeafNodeType &leaf, size_t leafIndex) const
Definition: VolumeToMesh.h:1124
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:757
int getValueDepth(const Coord &xyz) const
Definition: ValueAccessor.h:229
SeamWeights(const TreeT &distTree, const Int16TreeT &refSignTree, IntTreeT &refIdxTree, QuantizedPointList &points, double iso)
Definition: VolumeToMesh.h:1936
void operator()(LeafNodeType &signLeaf, size_t leafIndex) const
Definition: VolumeToMesh.h:1950
tree::ValueAccessor< IntTreeT > IntAccessorT
Definition: VolumeToMesh.h:901
PartGen(const LeafManagerT &leafs, size_t partitions, size_t activePart)
Definition: VolumeToMesh.h:2816
PolygonPoolList & polygonPoolList()
Definition: VolumeToMesh.h:3993
Definition: VolumeToMesh.h:3077
void join(PartGen &rhs)
Definition: VolumeToMesh.h:2807
tree::ValueAccessor< const Int16TreeT > Int16AccessorT
Definition: VolumeToMesh.h:1913
AdaptivePrimBuilder()
Definition: VolumeToMesh.h:2370
UniformPrimBuilder()
Definition: VolumeToMesh.h:2332
tree::LeafManager< const SrcTreeT > LeafManagerT
Definition: VolumeToMesh.h:2789
TreeT::template ValueConverter< float >::Type FloatTreeT
Definition: VolumeToMesh.h:2031
BoolTreeT & tree()
Definition: VolumeToMesh.h:3474
const char & triangleFlags(size_t n) const
Definition: VolumeToMesh.h:145
void operator()(LeafNodeType &leaf, size_t) const
Definition: VolumeToMesh.h:2998
PointListCopy(const PointList &pointsIn, std::vector< Vec3s > &pointsOut)
Definition: VolumeToMesh.h:3757
void setSpatialAdaptivity(const math::Transform &distGridXForm, const FloatGridT &adaptivityField)
Definition: VolumeToMesh.h:2106
Definition: VolumeToMesh.h:1634
void run(bool threaded=true)
Definition: VolumeToMesh.h:1708
tree::ValueAccessor< const IntTreeT > IntCAccessorT
Definition: VolumeToMesh.h:1641
void run(bool threaded=true)
Definition: VolumeToMesh.h:3090
void run(bool threaded=true)
Definition: VolumeToMesh.h:2917
tree::ValueAccessor< const IntTreeT > IntAccessorT
Definition: VolumeToMesh.h:1910
TreeT::ValueType ValueT
Definition: VolumeToMesh.h:2020
Definition: VolumeToMesh.h:3130
tree::ValueAccessor< const IntTreeT > IntAccessorT
Definition: VolumeToMesh.h:2621
OPENVDB_API const Coord COORD_OFFSETS[26]
coordinate offset table for neighboring voxels
CountPoints(std::vector< size_t > &pointList)
Definition: VolumeToMesh.h:1054
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:67
GenPoints(const LeafManagerT &signLeafs, const TreeT &distTree, IntTreeT &idxTree, PointList &points, std::vector< size_t > &indices, const math::Transform &xform, double iso)
Definition: VolumeToMesh.h:1687
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: VolumeToMesh.h:3101
BoolTreeT & tree()
Definition: VolumeToMesh.h:2800
void setRefData(const Int16TreeT *refSignTree=NULL, const TreeT *refDistTree=NULL, IntTreeT *refIdxTree=NULL, const QuantizedPointList *seamPoints=NULL, std::vector< unsigned char > *mSeamPointMaskPt=NULL)
Definition: VolumeToMesh.h:1717
Definition: VolumeToMesh.h:2991
TreeT::template ValueConverter< int >::Type IntTreeT
Definition: VolumeToMesh.h:900
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:109
tree::ValueAccessor< const SrcTreeT > SrcAccessorT
Definition: VolumeToMesh.h:3185
void setValueOn(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:246
Definition: VolumeToMesh.h:104
Definition: VolumeToMesh.h:397
TreeT::template ValueConverter< Int16 >::Type Int16TreeT
Definition: VolumeToMesh.h:1643
void setSpatialAdaptivity(const GridBase::ConstPtr &grid)
Definition: VolumeToMesh.h:4036
OPENVDB_STATIC_SPECIALIZATION void dilateVoxels(TreeType &tree, int count=1)
Definition: Morphology.h:364
Counts the total number of points per collapsed region.
Definition: VolumeToMesh.h:1111
PointList & pointList()
Definition: VolumeToMesh.h:3979
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:287
Definition: VolumeToMesh.h:397
tree::ValueAccessor< const BoolTreeT > BoolAccessorT
Definition: VolumeToMesh.h:2993
boost::scoped_array< uint32_t > QuantizedPointList
Definition: VolumeToMesh.h:1646
SrcTreeT::template ValueConverter< int >::Type IntTreeT
Definition: VolumeToMesh.h:3315
BoolTreeT & mask()
Definition: VolumeToMesh.h:2880
Definition: VolumeToMesh.h:1904
MergeVoxelRegions(const LeafManagerT &signLeafs, const Int16TreeT &signTree, const TreeT &distTree, IntTreeT &idxTree, ValueT iso, ValueT adaptivity)
Definition: VolumeToMesh.h:2075
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:347
SrcTreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: VolumeToMesh.h:3316
void constructPolygons(Int16 flags, Int16 refFlags, const Vec4i &offsets, const Coord &ijk, const SignAccT &signAcc, const IdxAccT &idxAcc, PrimBuilder &mesher, Index32 pointListSize)
Definition: VolumeToMesh.h:2465
tree::ValueAccessor< const TreeT > AccessorT
Definition: VolumeToMesh.h:2962
TreeT::template ValueConverter< Int16 >::Type Int16TreeT
Definition: VolumeToMesh.h:903
Definition: Mat4.h:51
Definition: VolumeToMesh.h:396
void clearTriangles()
Definition: VolumeToMesh.h:3885
void run(bool threaded=true)
Definition: VolumeToMesh.h:2840
Definition: VolumeToMesh.h:3461
Calculate an axis-aligned bounding box in index space from a bounding sphere in world space...
Definition: Transform.h:65
void operator()(const tbb::blocked_range< size_t > &) const
Definition: VolumeToMesh.h:1730
void operator()(LeafNodeType &leaf, size_t leafIndex) const
Definition: VolumeToMesh.h:1088
Definition: VolumeToMesh.h:2959
void run(bool threaded=true)
Definition: VolumeToMesh.h:2666