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