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