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