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