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