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