GCC Code Coverage Report

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