10 #ifndef OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
11 #define OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
18 #include <tbb/blocked_range.h>
19 #include <tbb/parallel_for.h>
20 #include <tbb/parallel_reduce.h>
21 #include <tbb/task_scheduler_init.h>
27 #include <type_traits>
51 template<
typename Gr
idType>
55 std::vector<Vec3s>& points,
56 std::vector<Vec4I>& quads,
72 template<
typename Gr
idType>
76 std::vector<Vec3s>& points,
77 std::vector<Vec3I>& triangles,
78 std::vector<Vec4I>& quads,
80 double adaptivity = 0.0,
97 inline PolygonPool(
const size_t numQuads,
const size_t numTriangles);
101 inline void resetQuads(
size_t size);
102 inline void clearQuads();
104 inline void resetTriangles(
size_t size);
105 inline void clearTriangles();
110 const size_t&
numQuads()
const {
return mNumQuads; }
125 const char&
quadFlags(
size_t n)
const {
return mQuadFlags[n]; }
134 inline bool trimQuads(
const size_t n,
bool reallocate =
false);
135 inline bool trimTrinagles(
const size_t n,
bool reallocate =
false);
141 size_t mNumQuads, mNumTriangles;
142 std::unique_ptr<openvdb::Vec4I[]> mQuads;
143 std::unique_ptr<openvdb::Vec3I[]> mTriangles;
144 std::unique_ptr<char[]> mQuadFlags, mTriangleFlags;
182 const std::vector<uint8_t>&
pointFlags()
const {
return mPointFlags; }
191 template<
typename InputGr
idType>
192 void operator()(
const InputGridType&);
245 size_t mPointListSize, mSeamPointListSize, mPolygonPoolListSize;
246 double mIsovalue, mPrimAdaptivity, mSecAdaptivity;
253 bool mInvertSurfaceMask, mRelaxDisorientedTriangles;
255 std::unique_ptr<uint32_t[]> mQuantizedSeamPoints;
256 std::vector<uint8_t> mPointFlags;
270 const std::vector<Vec3d>& points,
271 const std::vector<Vec3d>& normals)
277 if (points.empty())
return avgPos;
279 for (
size_t n = 0, N = points.size(); n < N; ++n) {
283 avgPos /= double(points.size());
287 double m00=0,m01=0,m02=0,
294 for (
size_t n = 0, N = points.size(); n < N; ++n) {
296 const Vec3d& n_ref = normals[n];
299 m00 += n_ref[0] * n_ref[0];
300 m11 += n_ref[1] * n_ref[1];
301 m22 += n_ref[2] * n_ref[2];
303 m01 += n_ref[0] * n_ref[1];
304 m02 += n_ref[0] * n_ref[2];
305 m12 += n_ref[1] * n_ref[2];
308 rhs += n_ref * n_ref.
dot(points[n] - avgPos);
333 Mat3d D = Mat3d::identity();
336 double tolerance =
std::max(std::abs(eigenValues[0]), std::abs(eigenValues[1]));
337 tolerance =
std::max(tolerance, std::abs(eigenValues[2]));
341 for (
int i = 0; i < 3; ++i ) {
342 if (std::abs(eigenValues[i]) < tolerance) {
346 D[i][i] = 1.0 / eigenValues[i];
353 return avgPos + pseudoInv * rhs;
367 namespace volume_to_mesh_internal {
369 template<
typename ValueType>
372 FillArray(ValueType* array,
const ValueType& v) : mArray(array), mValue(v) { }
374 void operator()(
const tbb::blocked_range<size_t>& range)
const {
375 const ValueType v = mValue;
376 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
386 template<
typename ValueType>
388 fillArray(ValueType* array,
const ValueType& val,
const size_t length)
390 const auto grainSize = std::max<size_t>(
391 length / tbb::task_scheduler_init::default_num_threads(), 1024);
392 const tbb::blocked_range<size_t> range(0, length, grainSize);
404 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,
405 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,
406 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,
407 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,
408 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,
409 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,
410 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,
411 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};
416 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,
417 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,
418 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,
419 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,
420 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,
421 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,
422 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,
423 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};
430 {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},
431 {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},
432 {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},
433 {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},
434 {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},
435 {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},
436 {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},
437 {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},
438 {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},
439 {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},
440 {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},
441 {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},
442 {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},
443 {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},
444 {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},
445 {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},
446 {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},
447 {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},
448 {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},
449 {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},
450 {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},
451 {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},
452 {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},
453 {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},
454 {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},
455 {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},
456 {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},
457 {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},
458 {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},
459 {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},
460 {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},
461 {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},
462 {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},
463 {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},
464 {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},
465 {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},
466 {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},
467 {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},
468 {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},
469 {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},
470 {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},
471 {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},
472 {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},
473 {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},
474 {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},
475 {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},
476 {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},
477 {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},
478 {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},
479 {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},
480 {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},
481 {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},
482 {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},
483 {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},
484 {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},
485 {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},
486 {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},
487 {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},
488 {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},
489 {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},
490 {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},
491 {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},
492 {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},
493 {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},
494 {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},
495 {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},
496 {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},
497 {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},
498 {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},
499 {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},
500 {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},
501 {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},
502 {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},
503 {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},
504 {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},
505 {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},
506 {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},
507 {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},
508 {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},
509 {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},
510 {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},
511 {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},
512 {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},
513 {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},
514 {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},
515 {0,0,0,0,0,0,0,0,0,0,0,0,0}};
524 double epsilon = 0.001)
527 Vec3d normal = (p2-p0).cross(p1-p3);
529 const Vec3d centroid = (p0 + p1 + p2 + p3);
530 const double d = centroid.
dot(normal) * 0.25;
534 double absDist = std::abs(p0.dot(normal) - d);
535 if (absDist > epsilon)
return false;
537 absDist = std::abs(p1.dot(normal) - d);
538 if (absDist > epsilon)
return false;
540 absDist = std::abs(p2.dot(normal) - d);
541 if (absDist > epsilon)
return false;
543 absDist = std::abs(p3.dot(normal) - d);
544 if (absDist > epsilon)
return false;
568 assert(!(v.x() > 1.0) && !(v.y() > 1.0) && !(v.z() > 1.0));
569 assert(!(v.x() < 0.0) && !(v.y() < 0.0) && !(v.z() < 0.0));
610 template<
typename AccessorT>
615 values[0] = accessor.getValue(ijk);
617 values[1] = accessor.getValue(ijk);
619 values[2] = accessor.getValue(ijk);
621 values[3] = accessor.getValue(ijk);
623 values[4] = accessor.getValue(ijk);
625 values[5] = accessor.getValue(ijk);
627 values[6] = accessor.getValue(ijk);
629 values[7] = accessor.getValue(ijk);
633 template<
typename LeafT>
638 values[0] = leaf.getValue(offset);
639 values[3] = leaf.getValue(offset + 1);
640 values[4] = leaf.getValue(offset + LeafT::DIM);
641 values[7] = leaf.getValue(offset + LeafT::DIM + 1);
642 values[1] = leaf.getValue(offset + (LeafT::DIM * LeafT::DIM));
643 values[2] = leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + 1);
644 values[5] = leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM);
645 values[6] = leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM + 1);
649 template<
typename ValueType>
662 return uint8_t(signs);
668 template<
typename AccessorT>
674 if (
isInsideValue(accessor.getValue(coord), iso)) signs |= 1u;
676 if (
isInsideValue(accessor.getValue(coord), iso)) signs |= 2u;
678 if (
isInsideValue(accessor.getValue(coord), iso)) signs |= 4u;
680 if (
isInsideValue(accessor.getValue(coord), iso)) signs |= 8u;
681 coord[1] += 1; coord[2] = ijk[2];
682 if (
isInsideValue(accessor.getValue(coord), iso)) signs |= 16u;
684 if (
isInsideValue(accessor.getValue(coord), iso)) signs |= 32u;
686 if (
isInsideValue(accessor.getValue(coord), iso)) signs |= 64u;
688 if (
isInsideValue(accessor.getValue(coord), iso)) signs |= 128u;
689 return uint8_t(signs);
695 template<
typename LeafT>
705 if (
isInsideValue(leaf.getValue(offset + 1), iso)) signs |= 8u;
708 if (
isInsideValue(leaf.getValue(offset + LeafT::DIM), iso)) signs |= 16u;
711 if (
isInsideValue(leaf.getValue(offset + LeafT::DIM + 1), iso)) signs |= 128u;
714 if (
isInsideValue(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) ), iso)) signs |= 2u;
717 if (
isInsideValue(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + 1), iso)) signs |= 4u;
720 if (
isInsideValue(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM), iso)) signs |= 32u;
723 if (
isInsideValue(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM + 1), iso)) signs |= 64u;
725 return uint8_t(signs);
731 template<
class AccessorT>
734 const AccessorT& acc,
Coord ijk,
typename AccessorT::ValueType iso)
767 template<
class AccessorT>
770 typename AccessorT::ValueType
isovalue,
const int dim)
783 coord[1] += dim; coord[2] = ijk[2];
794 if (p[0]) signs |= 1u;
795 if (p[1]) signs |= 2u;
796 if (p[2]) signs |= 4u;
797 if (p[3]) signs |= 8u;
798 if (p[4]) signs |= 16u;
799 if (p[5]) signs |= 32u;
800 if (p[6]) signs |= 64u;
801 if (p[7]) signs |= 128u;
807 int i = ijk[0], ip = ijk[0] + hDim, ipp = ijk[0] + dim;
808 int j = ijk[1], jp = ijk[1] + hDim, jpp = ijk[1] + dim;
809 int k = ijk[2], kp = ijk[2] + hDim, kpp = ijk[2] + dim;
812 coord.
reset(ip, j, k);
814 if (p[0] != m && p[1] != m)
return true;
817 coord.
reset(ipp, j, kp);
819 if (p[1] != m && p[2] != m)
return true;
822 coord.
reset(ip, j, kpp);
824 if (p[2] != m && p[3] != m)
return true;
827 coord.
reset(i, j, kp);
829 if (p[0] != m && p[3] != m)
return true;
832 coord.
reset(ip, jpp, k);
834 if (p[4] != m && p[5] != m)
return true;
837 coord.
reset(ipp, jpp, kp);
839 if (p[5] != m && p[6] != m)
return true;
842 coord.
reset(ip, jpp, kpp);
844 if (p[6] != m && p[7] != m)
return true;
847 coord.
reset(i, jpp, kp);
849 if (p[7] != m && p[4] != m)
return true;
852 coord.
reset(i, jp, k);
854 if (p[0] != m && p[4] != m)
return true;
857 coord.
reset(ipp, jp, k);
859 if (p[1] != m && p[5] != m)
return true;
862 coord.
reset(ipp, jp, kpp);
864 if (p[2] != m && p[6] != m)
return true;
868 coord.
reset(i, jp, kpp);
870 if (p[3] != m && p[7] != m)
return true;
876 coord.
reset(ip, jp, k);
878 if (p[0] != m && p[1] != m && p[4] != m && p[5] != m)
return true;
881 coord.
reset(ipp, jp, kp);
883 if (p[1] != m && p[2] != m && p[5] != m && p[6] != m)
return true;
886 coord.
reset(ip, jp, kpp);
888 if (p[2] != m && p[3] != m && p[6] != m && p[7] != m)
return true;
891 coord.
reset(i, jp, kp);
893 if (p[0] != m && p[3] != m && p[4] != m && p[7] != m)
return true;
896 coord.
reset(ip, j, kp);
898 if (p[0] != m && p[1] != m && p[2] != m && p[3] != m)
return true;
901 coord.
reset(ip, jpp, kp);
903 if (p[4] != m && p[5] != m && p[6] != m && p[7] != m)
return true;
906 coord.
reset(ip, jp, kp);
908 if (p[0] != m && p[1] != m && p[2] != m && p[3] != m &&
909 p[4] != m && p[5] != m && p[6] != m && p[7] != m)
return true;
918 template <
class LeafType>
922 Coord ijk, end = start;
927 for (ijk[0] = start[0]; ijk[0] < end[0]; ++ijk[0]) {
928 for (ijk[1] = start[1]; ijk[1] < end[1]; ++ijk[1]) {
929 for (ijk[2] = start[2]; ijk[2] < end[2]; ++ijk[2]) {
930 leaf.setValueOnly(ijk, regionId);
939 template <
class LeafType>
942 typename LeafType::ValueType::value_type adaptivity)
944 if (adaptivity < 1e-6)
return false;
946 using VecT =
typename LeafType::ValueType;
947 Coord ijk, end = start;
952 std::vector<VecT> norms;
953 for (ijk[0] = start[0]; ijk[0] < end[0]; ++ijk[0]) {
954 for (ijk[1] = start[1]; ijk[1] < end[1]; ++ijk[1]) {
955 for (ijk[2] = start[2]; ijk[2] < end[2]; ++ijk[2]) {
957 if(!leaf.isValueOn(ijk))
continue;
958 norms.push_back(leaf.getValue(ijk));
963 size_t N = norms.size();
964 for (
size_t ni = 0; ni < N; ++ni) {
965 VecT n_i = norms[ni];
966 for (
size_t nj = 0; nj < N; ++nj) {
967 VecT n_j = norms[nj];
968 if ((1.0 - n_i.dot(n_j)) > adaptivity)
return false;
979 inline double evalZeroCrossing(
double v0,
double v1,
double iso) {
return (iso - v0) / (v1 - v0); }
983 template<
typename LeafT>
987 values[0] = double(leaf.getValue(offset));
988 values[3] = double(leaf.getValue(offset + 1));
989 values[4] = double(leaf.getValue(offset + LeafT::DIM));
990 values[7] = double(leaf.getValue(offset + LeafT::DIM + 1));
991 values[1] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM)));
992 values[2] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + 1));
993 values[5] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM));
994 values[6] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM + 1));
999 template<
typename AccessorT>
1004 values[0] = double(acc.getValue(coord));
1007 values[1] = double(acc.getValue(coord));
1010 values[2] = double(acc.getValue(coord));
1013 values[3] = double(acc.getValue(coord));
1015 coord[1] += 1; coord[2] = ijk[2];
1016 values[4] = double(acc.getValue(coord));
1019 values[5] = double(acc.getValue(coord));
1022 values[6] = double(acc.getValue(coord));
1025 values[7] = double(acc.getValue(coord));
1032 unsigned char edgeGroup,
double iso)
1034 Vec3d avg(0.0, 0.0, 0.0);
1110 double w = 1.0 / double(samples);
1124 unsigned char signsMask,
unsigned char edgeGroup,
double iso)
1126 avg =
Vec3d(0.0, 0.0, 0.0);
1214 double w = 1.0 / double(samples);
1228 unsigned char signs,
unsigned char edgeGroup,
double iso)
1230 std::vector<Vec3d> samples;
1233 std::vector<double> weights;
1236 Vec3d avg(0.0, 0.0, 0.0);
1243 samples.push_back(avg);
1244 weights.push_back((avg-p).lengthSqr());
1252 samples.push_back(avg);
1253 weights.push_back((avg-p).lengthSqr());
1261 samples.push_back(avg);
1262 weights.push_back((avg-p).lengthSqr());
1270 samples.push_back(avg);
1271 weights.push_back((avg-p).lengthSqr());
1279 samples.push_back(avg);
1280 weights.push_back((avg-p).lengthSqr());
1288 samples.push_back(avg);
1289 weights.push_back((avg-p).lengthSqr());
1297 samples.push_back(avg);
1298 weights.push_back((avg-p).lengthSqr());
1306 samples.push_back(avg);
1307 weights.push_back((avg-p).lengthSqr());
1315 samples.push_back(avg);
1316 weights.push_back((avg-p).lengthSqr());
1324 samples.push_back(avg);
1325 weights.push_back((avg-p).lengthSqr());
1333 samples.push_back(avg);
1334 weights.push_back((avg-p).lengthSqr());
1342 samples.push_back(avg);
1343 weights.push_back((avg-p).lengthSqr());
1350 for (
size_t i = 0, I = weights.size(); i < I; ++i) {
1351 minWeight =
std::min(minWeight, weights[i]);
1352 maxWeight =
std::max(maxWeight, weights[i]);
1355 const double offset = maxWeight + minWeight * 0.1;
1356 for (
size_t i = 0, I = weights.size(); i < I; ++i) {
1357 weights[i] = offset - weights[i];
1361 double weightSum = 0.0;
1362 for (
size_t i = 0, I = weights.size(); i < I; ++i) {
1363 weightSum += weights[i];
1370 if (samples.size() > 1) {
1371 for (
size_t i = 0, I = samples.size(); i < I; ++i) {
1372 avg += samples[i] * (weights[i] / weightSum);
1375 avg = samples.front();
1386 const std::vector<double>& values,
unsigned char signs,
double iso)
1389 points.push_back(
computePoint(values, signs, uint8_t(n), iso));
1398 matchEdgeGroup(
unsigned char groupId,
unsigned char lhsSigns,
unsigned char rhsSigns)
1401 for (
size_t i = 1; i <= 12; ++i) {
1417 const std::vector<double>& lhsValues,
const std::vector<double>& rhsValues,
1418 unsigned char lhsSigns,
unsigned char rhsSigns,
1419 double iso,
size_t pointIdx,
const uint32_t * seamPointArray)
1421 for (
size_t n = 1, N =
sEdgeGroupTable[lhsSigns][0] + 1; n < N; ++n) {
1427 const unsigned char e = uint8_t(
id);
1428 const uint32_t& quantizedPoint = seamPointArray[pointIdx + (
id - 1)];
1433 weightedPointMask.push_back(
true);
1435 points.push_back(
computePoint(rhsValues, rhsSigns, e, iso));
1436 weightedPointMask.push_back(
false);
1440 points.push_back(
computePoint(lhsValues, lhsSigns, uint8_t(n), iso));
1441 weightedPointMask.push_back(
false);
1447 template <
typename InputTreeType>
1453 using Int16TreeType =
typename InputTreeType::template ValueConverter<Int16>::Type;
1460 const InputTreeType& inputTree,
1461 const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
1462 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
1463 const std::unique_ptr<
Index32[]>& leafNodeOffsets,
1467 void setRefData(
const InputTreeType& refInputTree,
1470 const uint32_t * quantizedSeamLinePoints,
1471 uint8_t * seamLinePointsFlags);
1473 void operator()(
const tbb::blocked_range<size_t>&)
const;
1476 Vec3s *
const mPoints;
1477 InputTreeType
const *
const mInputTree;
1480 Index32 const *
const mNodeOffsets;
1482 double const mIsovalue;
1484 InputTreeType
const * mRefInputTree;
1487 uint32_t
const * mQuantizedSeamLinePoints;
1488 uint8_t * mSeamLinePointsFlags;
1492 template <
typename InputTreeType>
1495 const InputTreeType& inputTree,
1496 const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
1497 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
1498 const std::unique_ptr<
Index32[]>& leafNodeOffsets,
1501 : mPoints(pointArray)
1502 , mInputTree(&inputTree)
1503 , mPointIndexNodes(pointIndexLeafNodes.empty() ? nullptr : &pointIndexLeafNodes.front())
1504 , mSignFlagsNodes(signFlagsLeafNodes.empty() ? nullptr : &signFlagsLeafNodes.front())
1505 , mNodeOffsets(leafNodeOffsets.get())
1508 , mRefInputTree(nullptr)
1509 , mRefPointIndexTree(nullptr)
1510 , mRefSignFlagsTree(nullptr)
1511 , mQuantizedSeamLinePoints(nullptr)
1512 , mSeamLinePointsFlags(nullptr)
1516 template <
typename InputTreeType>
1519 const InputTreeType& refInputTree,
1522 const uint32_t * quantizedSeamLinePoints,
1523 uint8_t * seamLinePointsFlags)
1525 mRefInputTree = &refInputTree;
1526 mRefPointIndexTree = &refPointIndexTree;
1527 mRefSignFlagsTree = &refSignFlagsTree;
1528 mQuantizedSeamLinePoints = quantizedSeamLinePoints;
1529 mSeamLinePointsFlags = seamLinePointsFlags;
1532 template <
typename InputTreeType>
1540 using IndexType =
typename Index32TreeType::ValueType;
1543 using IndexArrayMap = std::map<IndexType, IndexArray>;
1545 InputTreeAccessor inputAcc(*mInputTree);
1549 std::vector<Vec3d> points(4);
1550 std::vector<bool> weightedPointMask(4);
1551 std::vector<double> values(8), refValues(8);
1552 const double iso = mIsovalue;
1556 std::unique_ptr<InputTreeAccessor> refInputAcc;
1557 std::unique_ptr<Index32TreeAccessor> refPointIndexAcc;
1558 std::unique_ptr<Int16TreeAccessor> refSignFlagsAcc;
1560 const bool hasReferenceData = mRefInputTree && mRefPointIndexTree && mRefSignFlagsTree;
1562 if (hasReferenceData) {
1563 refInputAcc.reset(
new InputTreeAccessor(*mRefInputTree));
1564 refPointIndexAcc.reset(
new Index32TreeAccessor(*mRefPointIndexTree));
1565 refSignFlagsAcc.reset(
new Int16TreeAccessor(*mRefSignFlagsTree));
1568 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
1571 const Coord& origin = pointIndexNode.origin();
1581 if (hasReferenceData) {
1582 refInputNode = refInputAcc->probeConstLeaf(origin);
1583 refPointIndexNode = refPointIndexAcc->probeConstLeaf(origin);
1584 refSignFlagsNode = refSignFlagsAcc->probeConstLeaf(origin);
1587 IndexType pointOffset = IndexType(mNodeOffsets[n]);
1588 IndexArrayMap regions;
1590 for (
auto it = pointIndexNode.beginValueOn(); it; ++it) {
1591 const Index offset = it.pos();
1593 const IndexType
id = it.getValue();
1596 regions[id].push_back(offset);
1601 pointIndexNode.setValueOnly(offset, pointOffset);
1603 const Int16 flags = signFlagsNode.getValue(offset);
1604 uint8_t signs = uint8_t(
SIGNS & flags);
1605 uint8_t refSigns = 0;
1607 if ((flags &
SEAM) && refPointIndexNode && refSignFlagsNode) {
1608 if (refSignFlagsNode->isValueOn(offset)) {
1609 refSigns = uint8_t(
SIGNS & refSignFlagsNode->getValue(offset));
1613 ijk = Index32LeafNodeType::offsetToLocalCoord(offset);
1615 const bool inclusiveCell = inputNode &&
1616 ijk[0] < int(Index32LeafNodeType::DIM - 1) &&
1617 ijk[1] < int(Index32LeafNodeType::DIM - 1) &&
1618 ijk[2] < int(Index32LeafNodeType::DIM - 1);
1626 weightedPointMask.clear();
1628 if (refSigns == 0) {
1633 if (inclusiveCell && refInputNode) {
1638 computeCellPoints(points, weightedPointMask, values, refValues, signs, refSigns,
1639 iso, refPointIndexNode->getValue(offset), mQuantizedSeamLinePoints);
1642 xyz[0] = double(ijk[0]);
1643 xyz[1] = double(ijk[1]);
1644 xyz[2] = double(ijk[2]);
1646 for (
size_t i = 0, I = points.size(); i < I; ++i) {
1648 Vec3d& point = points[i];
1651 if (!std::isfinite(point[0]) ||
1652 !std::isfinite(point[1]) ||
1653 !std::isfinite(point[2]))
1656 "VolumeToMesh encountered NaNs or infs in the input VDB!"
1657 " Hint: Check the input and consider using the \"Diagnostics\" tool "
1658 "to detect and resolve the NaNs.");
1662 point = mTransform.indexToWorld(point);
1664 Vec3s& pos = mPoints[pointOffset];
1665 pos[0] = float(point[0]);
1666 pos[1] = float(point[1]);
1667 pos[2] = float(point[2]);
1669 if (mSeamLinePointsFlags && !weightedPointMask.empty() && weightedPointMask[i]) {
1670 mSeamLinePointsFlags[pointOffset] = uint8_t(1);
1678 for (
typename IndexArrayMap::iterator it = regions.begin(); it != regions.end(); ++it) {
1680 Vec3d avg(0.0), point(0.0);
1684 for (
size_t i = 0, I = voxels.size(); i < I; ++i) {
1686 const Index offset = voxels[i];
1687 ijk = Index32LeafNodeType::offsetToLocalCoord(offset);
1689 const bool inclusiveCell = inputNode &&
1690 ijk[0] < int(Index32LeafNodeType::DIM - 1) &&
1691 ijk[1] < int(Index32LeafNodeType::DIM - 1) &&
1692 ijk[2] < int(Index32LeafNodeType::DIM - 1);
1696 pointIndexNode.setValueOnly(offset, pointOffset);
1698 uint8_t signs = uint8_t(
SIGNS & signFlagsNode.getValue(offset));
1706 avg[0] += double(ijk[0]) + points[0][0];
1707 avg[1] += double(ijk[1]) + points[0][1];
1708 avg[2] += double(ijk[2]) + points[0][2];
1714 double w = 1.0 / double(count);
1720 avg = mTransform.indexToWorld(avg);
1722 Vec3s& pos = mPoints[pointOffset];
1723 pos[0] = float(avg[0]);
1724 pos[1] = float(avg[1]);
1725 pos[2] = float(avg[2]);
1736 template <
typename InputTreeType>
1742 using Int16TreeType =
typename InputTreeType::template ValueConverter<Int16>::Type;
1749 const InputTreeType& inputTree,
1752 uint32_t * quantizedPoints,
1754 : mSignFlagsNodes(signFlagsLeafNodes.empty() ? nullptr : &signFlagsLeafNodes.front())
1755 , mInputTree(&inputTree)
1756 , mRefPointIndexTree(&refPointIndexTree)
1757 , mRefSignFlagsTree(&refSignFlagsTree)
1758 , mQuantizedPoints(quantizedPoints)
1769 std::vector<double> values(8);
1770 const double iso = double(mIsovalue);
1774 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
1777 const Coord& origin = signFlagsNode.origin();
1780 if (!refSignNode)
continue;
1784 if (!refPointIndexNode)
continue;
1788 for (
typename Int16LeafNodeType::ValueOnCIter it = signFlagsNode.cbeginValueOn();
1791 const Index offset = it.pos();
1793 ijk = Index32LeafNodeType::offsetToLocalCoord(offset);
1795 const bool inclusiveCell = inputNode &&
1796 ijk[0] < int(Index32LeafNodeType::DIM - 1) &&
1797 ijk[1] < int(Index32LeafNodeType::DIM - 1) &&
1798 ijk[2] < int(Index32LeafNodeType::DIM - 1);
1802 if ((it.getValue() &
SEAM) && refSignNode->isValueOn(offset)) {
1804 uint8_t lhsSigns = uint8_t(
SIGNS & it.getValue());
1805 uint8_t rhsSigns = uint8_t(
SIGNS & refSignNode->getValue(offset));
1808 if (inclusiveCell) {
1815 for (
unsigned i = 1, I =
sEdgeGroupTable[lhsSigns][0] + 1; i < I; ++i) {
1821 uint32_t& data = mQuantizedPoints[
1822 refPointIndexNode->getValue(offset) + (
id - 1)];
1827 pos, values, lhsSigns, rhsSigns, uint8_t(i), iso);
1845 InputTreeType
const *
const mInputTree;
1848 uint32_t *
const mQuantizedPoints;
1853 template <
typename TreeType>
1859 const TreeType& refSignFlagsTree)
1860 : mSignFlagsNodes(signFlagsLeafNodes.empty() ? nullptr : &signFlagsLeafNodes.front())
1861 , mRefSignFlagsTree(&refSignFlagsTree)
1869 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
1872 const Coord& origin = signFlagsNode.origin();
1875 if (!refSignNode)
continue;
1877 for (
auto it = signFlagsNode.cbeginValueOn(); it; ++it) {
1878 const Index offset = it.pos();
1880 uint8_t rhsSigns = uint8_t(refSignNode->getValue(offset) &
SIGNS);
1884 const typename LeafNodeType::ValueType value = it.getValue();
1885 uint8_t lhsSigns = uint8_t(value &
SIGNS);
1887 if (rhsSigns != lhsSigns) {
1888 signFlagsNode.setValueOnly(offset, value |
SEAM);
1899 TreeType
const *
const mRefSignFlagsTree;
1903 template <
typename BoolTreeType,
typename SignDataType>
1912 const BoolTreeType& maskTree)
1913 : mSignFlagsNodes(signFlagsLeafNodes.empty() ? nullptr : &signFlagsLeafNodes.front())
1914 , mMaskTree(&maskTree)
1922 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
1925 const Coord& origin = signFlagsNode.origin();
1928 if (!maskNode)
continue;
1930 using ValueOnCIter =
typename SignDataLeafNodeType::ValueOnCIter;
1932 for (ValueOnCIter it = signFlagsNode.cbeginValueOn(); it; ++it) {
1933 const Index offset = it.pos();
1935 if (maskNode->isValueOn(offset)) {
1936 signFlagsNode.setValueOnly(offset, it.getValue() |
SEAM);
1946 BoolTreeType
const *
const mMaskTree;
1950 template <
typename TreeType>
1954 using BoolTreeType =
typename TreeType::template ValueConverter<bool>::Type;
1957 const TreeType& signFlagsTree,
1959 : mSignFlagsNodes(signFlagsLeafNodes.empty() ? nullptr : &signFlagsLeafNodes.front())
1960 , mSignFlagsTree(&signFlagsTree)
1967 : mSignFlagsNodes(rhs.mSignFlagsNodes)
1968 , mSignFlagsTree(rhs.mSignFlagsTree)
1978 using ValueOnCIter =
typename LeafNodeType::ValueOnCIter;
1979 using ValueType =
typename LeafNodeType::ValueType;
1985 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
1990 for (ValueOnCIter it = signFlagsNode.cbeginValueOn(); it; ++it) {
1992 const ValueType flags = it.getValue();
1994 if (!(flags &
SEAM) && (flags &
EDGES)) {
1996 ijk = it.getCoord();
1998 bool isSeamLineVoxel =
false;
2000 if (flags &
XEDGE) {
2004 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) &
SEAM);
2006 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) &
SEAM);
2010 if (!isSeamLineVoxel && flags &
YEDGE) {
2012 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) &
SEAM);
2014 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) &
SEAM);
2016 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) &
SEAM);
2020 if (!isSeamLineVoxel && flags &
ZEDGE) {
2022 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) &
SEAM);
2024 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) &
SEAM);
2026 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) &
SEAM);
2030 if (isSeamLineVoxel) {
2031 maskAcc.
setValue(it.getCoord(),
true);
2041 TreeType
const *
const mSignFlagsTree;
2047 template<
typename SignDataTreeType>
2051 using SignDataType =
typename SignDataTreeType::ValueType;
2052 using SignDataLeafNodeType =
typename SignDataTreeType::LeafNodeType;
2053 using BoolTreeType =
typename SignDataTreeType::template ValueConverter<bool>::Type;
2055 std::vector<SignDataLeafNodeType*> signFlagsLeafNodes;
2056 signFlagsTree.getNodes(signFlagsLeafNodes);
2058 const tbb::blocked_range<size_t> nodeRange(0, signFlagsLeafNodes.size());
2060 tbb::parallel_for(nodeRange,
2063 BoolTreeType seamLineMaskTree(
false);
2066 maskSeamLine(signFlagsLeafNodes, signFlagsTree, seamLineMaskTree);
2068 tbb::parallel_reduce(nodeRange, maskSeamLine);
2070 tbb::parallel_for(nodeRange,
2078 template <
typename InputGr
idType>
2085 using FloatTreeType =
typename InputTreeType::template ValueConverter<float>::Type;
2089 using Int16TreeType =
typename InputTreeType::template ValueConverter<Int16>::Type;
2095 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
2100 const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
2101 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
2104 bool invertSurfaceOrientation);
2108 mSpatialAdaptivityTree = &grid.
tree();
2109 mSpatialAdaptivityTransform = &grid.transform();
2119 mRefSignFlagsTree = &signFlagsData;
2120 mInternalAdaptivity = internalAdaptivity;
2123 void operator()(
const tbb::blocked_range<size_t>&)
const;
2134 float mSurfaceAdaptivity, mInternalAdaptivity;
2135 bool mInvertSurfaceOrientation;
2144 template <
typename InputGr
idType>
2146 const InputGridType& inputGrid,
2148 const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
2149 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
2152 bool invertSurfaceOrientation)
2153 : mInputTree(&inputGrid.tree())
2154 , mInputTransform(&inputGrid.transform())
2155 , mPointIndexTree(&pointIndexTree)
2156 , mPointIndexNodes(pointIndexLeafNodes.empty() ? nullptr : &pointIndexLeafNodes.front())
2157 , mSignFlagsNodes(signFlagsLeafNodes.empty() ? nullptr : &signFlagsLeafNodes.front())
2159 , mSurfaceAdaptivity(adaptivity)
2160 , mInternalAdaptivity(adaptivity)
2161 , mInvertSurfaceOrientation(invertSurfaceOrientation)
2162 , mSpatialAdaptivityTree(nullptr)
2163 , mMaskTree(nullptr)
2164 , mRefSignFlagsTree(nullptr)
2165 , mSpatialAdaptivityTransform(nullptr)
2170 template <
typename InputGr
idType>
2175 using Vec3sLeafNodeType =
typename InputLeafNodeType::template ValueConverter<Vec3sType>::Type;
2183 std::unique_ptr<FloatTreeAccessor> spatialAdaptivityAcc;
2184 if (mSpatialAdaptivityTree && mSpatialAdaptivityTransform) {
2185 spatialAdaptivityAcc.reset(
new FloatTreeAccessor(*mSpatialAdaptivityTree));
2188 std::unique_ptr<BoolTreeAccessor> maskAcc;
2190 maskAcc.reset(
new BoolTreeAccessor(*mMaskTree));
2193 std::unique_ptr<Int16TreeAccessor> refSignFlagsAcc;
2194 if (mRefSignFlagsTree) {
2195 refSignFlagsAcc.reset(
new Int16TreeAccessor(*mRefSignFlagsTree));
2198 InputTreeAccessor inputAcc(*mInputTree);
2199 Index32TreeAccessor pointIndexAcc(*mPointIndexTree);
2203 const bool invertGradientDir = mInvertSurfaceOrientation || isBoolValue<InputValueType>();
2204 std::unique_ptr<Vec3sLeafNodeType> gradientNode;
2207 const int LeafDim = InputLeafNodeType::DIM;
2209 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
2211 mask.setValuesOff();
2216 const Coord& origin = pointIndexNode.origin();
2218 end[0] = origin[0] + LeafDim;
2219 end[1] = origin[1] + LeafDim;
2220 end[2] = origin[2] + LeafDim;
2225 if (maskLeaf !=
nullptr) {
2226 for (
typename BoolLeafNodeType::ValueOnCIter it = maskLeaf->cbeginValueOn();
2229 mask.setActiveState(it.getCoord() & ~1u,
true);
2234 float adaptivity = (refSignFlagsAcc && !refSignFlagsAcc->probeConstLeaf(origin)) ?
2235 mInternalAdaptivity : mSurfaceAdaptivity;
2237 bool useGradients = adaptivity < 1.0f;
2242 if (spatialAdaptivityAcc) {
2243 useGradients =
false;
2244 for (
Index offset = 0; offset < FloatLeafNodeType::NUM_VALUES; ++offset) {
2245 ijk = adaptivityLeaf.offsetToGlobalCoord(offset);
2246 ijk = mSpatialAdaptivityTransform->worldToIndexCellCentered(
2247 mInputTransform->indexToWorld(ijk));
2248 float weight = spatialAdaptivityAcc->getValue(ijk);
2249 float adaptivityValue = weight * adaptivity;
2250 if (adaptivityValue < 1.0f) useGradients =
true;
2251 adaptivityLeaf.setValueOnly(offset, adaptivityValue);
2256 for (
auto it = signFlagsNode.cbeginValueOn(); it; ++it) {
2257 const Int16 flags = it.getValue();
2258 const unsigned char signs =
static_cast<unsigned char>(
SIGNS & int(flags));
2262 mask.setActiveState(it.getCoord() & ~1u,
true);
2264 }
else if (flags &
EDGES) {
2266 bool maskRegion =
false;
2268 ijk = it.getCoord();
2269 if (!pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2271 if (!maskRegion && flags &
XEDGE) {
2273 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2275 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2277 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2281 if (!maskRegion && flags &
YEDGE) {
2283 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2285 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2287 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2291 if (!maskRegion && flags &
ZEDGE) {
2293 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2295 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2297 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2302 mask.setActiveState(it.getCoord() & ~1u,
true);
2309 for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
2310 for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
2311 for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
2312 if (!mask.isValueOn(ijk) &&
isNonManifold(inputAcc, ijk, mIsovalue, dim)) {
2313 mask.setActiveState(ijk,
true);
2324 gradientNode->setValuesOff();
2326 gradientNode.
reset(
new Vec3sLeafNodeType());
2329 for (
auto it = signFlagsNode.cbeginValueOn(); it; ++it) {
2330 ijk = it.getCoord();
2331 if (!mask.isValueOn(ijk & ~1u)) {
2335 if (invertGradientDir) {
2339 gradientNode->setValueOn(it.pos(), dir);
2346 for ( ; dim <= LeafDim; dim = dim << 1) {
2347 const unsigned coordMask = ~((dim << 1) - 1);
2348 for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
2349 for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
2350 for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
2352 adaptivity = adaptivityLeaf.getValue(ijk);
2354 if (mask.isValueOn(ijk)
2356 || (useGradients && !
isMergable(*gradientNode, ijk, dim, adaptivity)))
2358 mask.setActiveState(ijk & coordMask,
true);
2360 mergeVoxels(pointIndexNode, ijk, dim, regionId++);
2380 mPolygonPool = &quadPool;
2385 template<
typename IndexType>
2389 mPolygonPool->
quad(mIdx) = verts;
2419 mPolygonPool = &polygonPool;
2427 template<
typename IndexType>
2430 if (verts[0] != verts[1] && verts[0] != verts[2] && verts[0] != verts[3]
2431 && verts[1] != verts[2] && verts[1] != verts[3] && verts[2] != verts[3]) {
2432 mPolygonPool->
quadFlags(mQuadIdx) = flags;
2433 addQuad(verts, reverse);
2435 verts[0] == verts[3] &&
2436 verts[1] != verts[2] &&
2437 verts[1] != verts[0] &&
2438 verts[2] != verts[0]) {
2440 addTriangle(verts[0], verts[1], verts[2], reverse);
2442 verts[1] == verts[2] &&
2443 verts[0] != verts[3] &&
2444 verts[0] != verts[1] &&
2445 verts[3] != verts[1]) {
2447 addTriangle(verts[0], verts[1], verts[3], reverse);
2449 verts[0] == verts[1] &&
2450 verts[2] != verts[3] &&
2451 verts[2] != verts[0] &&
2452 verts[3] != verts[0]) {
2454 addTriangle(verts[0], verts[2], verts[3], reverse);
2456 verts[2] == verts[3] &&
2457 verts[0] != verts[1] &&
2458 verts[0] != verts[2] &&
2459 verts[1] != verts[2]) {
2461 addTriangle(verts[0], verts[1], verts[2], reverse);
2468 mPolygonPool->
trimQuads(mQuadIdx,
true);
2474 template<
typename IndexType>
2478 mPolygonPool->
quad(mQuadIdx) = verts;
2480 Vec4I& quad = mPolygonPool->
quad(mQuadIdx);
2489 void addTriangle(
unsigned v0,
unsigned v1,
unsigned v2,
bool reverse)
2505 size_t mQuadIdx, mTriangleIdx;
2506 PolygonPool *mPolygonPool;
2510 template<
typename SignAccT,
typename IdxAccT,
typename PrimBuilder>
2513 bool invertSurfaceOrientation,
2516 const Vec3i& offsets,
2518 const SignAccT& signAcc,
2519 const IdxAccT& idxAcc,
2520 PrimBuilder& mesher)
2522 using IndexType =
typename IdxAccT::ValueType;
2525 const bool isActive = idxAcc.probeValue(ijk, v0);
2532 bool isInside = flags &
INSIDE;
2534 isInside = invertSurfaceOrientation ? !isInside : isInside;
2539 if (flags &
XEDGE) {
2541 quad[0] = v0 + offsets[0];
2545 bool activeValues = idxAcc.probeValue(coord, quad[1]);
2546 uint8_t cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2551 activeValues = activeValues && idxAcc.probeValue(coord, quad[2]);
2552 cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2557 activeValues = activeValues && idxAcc.probeValue(coord, quad[3]);
2558 cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2562 mesher.addPrim(quad, isInside, tag[
bool(refFlags &
XEDGE)]);
2569 if (flags &
YEDGE) {
2571 quad[0] = v0 + offsets[1];
2575 bool activeValues = idxAcc.probeValue(coord, quad[1]);
2576 uint8_t cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2581 activeValues = activeValues && idxAcc.probeValue(coord, quad[2]);
2582 cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2587 activeValues = activeValues && idxAcc.probeValue(coord, quad[3]);
2588 cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2592 mesher.addPrim(quad, isInside, tag[
bool(refFlags &
YEDGE)]);
2599 if (flags &
ZEDGE) {
2601 quad[0] = v0 + offsets[2];
2605 bool activeValues = idxAcc.probeValue(coord, quad[1]);
2606 uint8_t cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2611 activeValues = activeValues && idxAcc.probeValue(coord, quad[2]);
2612 cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2617 activeValues = activeValues && idxAcc.probeValue(coord, quad[3]);
2618 cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2622 mesher.addPrim(quad, !isInside, tag[
bool(refFlags &
ZEDGE)]);
2631 template<
typename InputTreeType>
2635 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
2640 : mInputTree(&inputTree)
2644 , mTileArray(tileArray)
2649 : mInputTree(rhs.mInputTree)
2650 , mIsovalue(rhs.mIsovalue)
2653 , mTileArray(rhs.mTileArray)
2659 void operator()(
const tbb::blocked_range<size_t>&);
2662 InputTreeType
const *
const mInputTree;
2666 Vec4i const *
const mTileArray;
2670 template<
typename InputTreeType>
2679 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
2681 const Vec4i& tile = mTileArray[n];
2683 bbox.
min()[0] = tile[0];
2684 bbox.
min()[1] = tile[1];
2685 bbox.
min()[2] = tile[2];
2700 bool processRegion =
true;
2705 if (processRegion) {
2708 region.
min()[0] = region.
max()[0] = ijk[0];
2709 mMask->fill(region,
false);
2716 processRegion =
true;
2718 processRegion = (!inputTreeAcc.
probeValue(ijk, value)
2722 if (processRegion) {
2725 region.
min()[0] = region.
max()[0] = ijk[0];
2726 mMask->fill(region,
false);
2736 processRegion =
true;
2741 if (processRegion) {
2744 region.
min()[1] = region.
max()[1] = ijk[1];
2745 mMask->fill(region,
false);
2752 processRegion =
true;
2754 processRegion = (!inputTreeAcc.
probeValue(ijk, value)
2758 if (processRegion) {
2761 region.
min()[1] = region.
max()[1] = ijk[1];
2762 mMask->fill(region,
false);
2772 processRegion =
true;
2777 if (processRegion) {
2780 region.
min()[2] = region.
max()[2] = ijk[2];
2781 mMask->fill(region,
false);
2787 processRegion =
true;
2789 processRegion = (!inputTreeAcc.
probeValue(ijk, value)
2793 if (processRegion) {
2796 region.
min()[2] = region.
max()[2] = ijk[2];
2797 mMask->fill(region,
false);
2803 template<
typename InputTreeType>
2806 typename InputTreeType::template ValueConverter<bool>::Type& mask)
2808 typename InputTreeType::ValueOnCIter tileIter(inputTree);
2809 tileIter.setMaxDepth(InputTreeType::ValueOnCIter::LEAF_DEPTH - 1);
2811 size_t tileCount = 0;
2812 for ( ; tileIter; ++tileIter) {
2816 if (tileCount > 0) {
2817 std::unique_ptr<Vec4i[]> tiles(
new Vec4i[tileCount]);
2822 tileIter = inputTree.cbeginValueOn();
2823 tileIter.setMaxDepth(InputTreeType::ValueOnCIter::LEAF_DEPTH - 1);
2825 for (; tileIter; ++tileIter) {
2826 Vec4i& tile = tiles[index++];
2827 tileIter.getBoundingBox(bbox);
2828 tile[0] = bbox.
min()[0];
2829 tile[1] = bbox.
min()[1];
2830 tile[2] = bbox.
min()[2];
2831 tile[3] = bbox.
max()[0] - bbox.
min()[0];
2835 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, tileCount), op);
2848 : mPointsIn(pointsIn) , mPointsOut(pointsOut)
2854 for (
size_t n = range.begin(); n < range.end(); ++n) {
2855 mPointsOut[n] = mPointsIn[n];
2861 std::vector<Vec3s>& mPointsOut;
2873 template<
typename LeafNodeType>
2912 IndexVector mCore, mMinX, mMaxX, mMinY, mMaxY, mMinZ, mMaxZ,
2913 mInternalNeighborsX, mInternalNeighborsY, mInternalNeighborsZ;
2917 template<
typename LeafNodeType>
2923 mCore.reserve((LeafNodeType::DIM - 2) * (LeafNodeType::DIM - 2));
2925 for (
Index x = 1; x < (LeafNodeType::DIM - 1); ++x) {
2926 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2927 for (
Index y = 1; y < (LeafNodeType::DIM - 1); ++y) {
2928 const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
2929 for (
Index z = 1; z < (LeafNodeType::DIM - 1); ++z) {
2930 mCore.push_back(offsetXY + z);
2936 mInternalNeighborsX.clear();
2937 mInternalNeighborsX.reserve(LeafNodeType::SIZE - (LeafNodeType::DIM * LeafNodeType::DIM));
2939 for (
Index x = 0; x < (LeafNodeType::DIM - 1); ++x) {
2940 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2941 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
2942 const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
2943 for (
Index z = 0; z < LeafNodeType::DIM; ++z) {
2944 mInternalNeighborsX.push_back(offsetXY + z);
2950 mInternalNeighborsY.clear();
2951 mInternalNeighborsY.reserve(LeafNodeType::SIZE - (LeafNodeType::DIM * LeafNodeType::DIM));
2953 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
2954 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2955 for (
Index y = 0; y < (LeafNodeType::DIM - 1); ++y) {
2956 const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
2957 for (
Index z = 0; z < LeafNodeType::DIM; ++z) {
2958 mInternalNeighborsY.push_back(offsetXY + z);
2964 mInternalNeighborsZ.clear();
2965 mInternalNeighborsZ.reserve(LeafNodeType::SIZE - (LeafNodeType::DIM * LeafNodeType::DIM));
2967 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
2968 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2969 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
2970 const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
2971 for (
Index z = 0; z < (LeafNodeType::DIM - 1); ++z) {
2972 mInternalNeighborsZ.push_back(offsetXY + z);
2979 mMinX.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
2981 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
2982 const Index offsetXY = (y << LeafNodeType::LOG2DIM);
2983 for (
Index z = 0; z < LeafNodeType::DIM; ++z) {
2984 mMinX.push_back(offsetXY + z);
2991 mMaxX.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
2993 const Index offsetX = (LeafNodeType::DIM - 1) << (2 * LeafNodeType::LOG2DIM);
2994 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
2995 const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
2996 for (
Index z = 0; z < LeafNodeType::DIM; ++z) {
2997 mMaxX.push_back(offsetXY + z);
3004 mMinY.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
3006 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
3007 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
3008 for (
Index z = 0; z < (LeafNodeType::DIM - 1); ++z) {
3009 mMinY.push_back(offsetX + z);
3016 mMaxY.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
3018 const Index offsetY = (LeafNodeType::DIM - 1) << LeafNodeType::LOG2DIM;
3019 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
3020 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
3021 for (
Index z = 0; z < (LeafNodeType::DIM - 1); ++z) {
3022 mMaxY.push_back(offsetX + offsetY + z);
3029 mMinZ.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
3031 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
3032 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
3033 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
3034 const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
3035 mMinZ.push_back(offsetXY);
3042 mMaxZ.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
3044 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
3045 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
3046 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
3047 const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
3048 mMaxZ.push_back(offsetXY + (LeafNodeType::DIM - 1));
3059 template<
typename AccessorT,
int _AXIS>
3069 acc.setActiveState(ijk);
3071 acc.setActiveState(ijk);
3073 acc.setActiveState(ijk);
3075 acc.setActiveState(ijk);
3076 }
else if (_AXIS == 1) {
3077 acc.setActiveState(ijk);
3079 acc.setActiveState(ijk);
3081 acc.setActiveState(ijk);
3083 acc.setActiveState(ijk);
3085 acc.setActiveState(ijk);
3087 acc.setActiveState(ijk);
3089 acc.setActiveState(ijk);
3091 acc.setActiveState(ijk);
3100 template<
typename VoxelEdgeAcc,
typename LeafNode>
3108 if (VoxelEdgeAcc::AXIS == 0) {
3109 nvo = LeafNode::DIM * LeafNode::DIM;
3111 }
else if (VoxelEdgeAcc::AXIS == 1) {
3112 nvo = LeafNode::DIM;
3116 for (
size_t n = 0, N = offsets->size(); n < N; ++n) {
3117 const Index& pos = (*offsets)[n];
3118 bool isActive = leafnode.isValueOn(pos) || leafnode.isValueOn(pos + nvo);
3119 if (isActive && (
isInsideValue(leafnode.getValue(pos), iso) !=
3121 edgeAcc.set(leafnode.offsetToGlobalCoord(pos));
3130 template<
typename LeafNode,
typename TreeAcc,
typename VoxelEdgeAcc>
3135 const std::vector<Index>* lhsOffsets = &voxels.
maxX();
3136 const std::vector<Index>* rhsOffsets = &voxels.
minX();
3137 Coord ijk = lhsNode.origin();
3139 if (VoxelEdgeAcc::AXIS == 0) {
3140 ijk[0] += LeafNode::DIM;
3141 }
else if (VoxelEdgeAcc::AXIS == 1) {
3142 ijk[1] += LeafNode::DIM;
3143 lhsOffsets = &voxels.
maxY();
3144 rhsOffsets = &voxels.
minY();
3145 }
else if (VoxelEdgeAcc::AXIS == 2) {
3146 ijk[2] += LeafNode::DIM;
3147 lhsOffsets = &voxels.
maxZ();
3148 rhsOffsets = &voxels.
minZ();
3151 typename LeafNode::ValueType value;
3152 const LeafNode* rhsNodePt = acc.probeConstLeaf(ijk);
3155 for (
size_t n = 0, N = lhsOffsets->size(); n < N; ++n) {
3156 const Index& pos = (*lhsOffsets)[n];
3157 bool isActive = lhsNode.isValueOn(pos) || rhsNodePt->isValueOn((*rhsOffsets)[n]);
3158 if (isActive && (
isInsideValue(lhsNode.getValue(pos), iso) !=
3159 isInsideValue(rhsNodePt->getValue((*rhsOffsets)[n]), iso))) {
3160 edgeAcc.set(lhsNode.offsetToGlobalCoord(pos));
3163 }
else if (!acc.probeValue(ijk, value)) {
3165 for (
size_t n = 0, N = lhsOffsets->size(); n < N; ++n) {
3166 const Index& pos = (*lhsOffsets)[n];
3167 if (lhsNode.isValueOn(pos) && (inside !=
isInsideValue(lhsNode.getValue(pos), iso))) {
3168 edgeAcc.set(lhsNode.offsetToGlobalCoord(pos));
3178 template<
typename LeafNode,
typename TreeAcc,
typename VoxelEdgeAcc>
3183 Coord ijk = leafnode.origin();
3184 if (VoxelEdgeAcc::AXIS == 0) --ijk[0];
3185 else if (VoxelEdgeAcc::AXIS == 1) --ijk[1];
3186 else if (VoxelEdgeAcc::AXIS == 2) --ijk[2];
3188 typename LeafNode::ValueType value;
3189 if (!acc.probeConstLeaf(ijk) && !acc.probeValue(ijk, value)) {
3196 for (
size_t n = 0, N = offsets->size(); n < N; ++n) {
3198 const Index& pos = (*offsets)[n];
3199 if (leafnode.isValueOn(pos)
3202 ijk = leafnode.offsetToGlobalCoord(pos);
3203 if (VoxelEdgeAcc::AXIS == 0) --ijk[0];
3204 else if (VoxelEdgeAcc::AXIS == 1) --ijk[1];
3205 else if (VoxelEdgeAcc::AXIS == 2) --ijk[2];
3215 template<
typename InputTreeType>
3221 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
3224 const InputTreeType& inputTree,
3225 const std::vector<const InputLeafNodeType*>& inputLeafNodes,
3230 void operator()(
const tbb::blocked_range<size_t>&);
3232 mIntersectionAccessor.
tree().merge(rhs.mIntersectionAccessor.
tree());
3249 template<
typename InputTreeType>
3251 const InputTreeType& inputTree,
3252 const std::vector<const InputLeafNodeType*>& inputLeafNodes,
3255 : mInputAccessor(inputTree)
3256 , mInputNodes(inputLeafNodes.empty() ? nullptr : &inputLeafNodes.front())
3257 , mIntersectionTree(false)
3258 , mIntersectionAccessor(intersectionTree)
3260 , mOffsets(&mOffsetData)
3267 template<
typename InputTreeType>
3270 : mInputAccessor(rhs.mInputAccessor.tree())
3271 , mInputNodes(rhs.mInputNodes)
3272 , mIntersectionTree(false)
3273 , mIntersectionAccessor(mIntersectionTree)
3275 , mOffsets(rhs.mOffsets)
3276 , mIsovalue(rhs.mIsovalue)
3281 template<
typename InputTreeType>
3289 for (
size_t n = range.begin(); n != range.end(); ++n) {
3320 template<
typename InputTreeType>
3323 typename InputTreeType::template ValueConverter<bool>::Type& intersectionTree,
3324 const InputTreeType& inputTree,
3325 typename InputTreeType::ValueType
isovalue)
3327 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
3329 std::vector<const InputLeafNodeType*> inputLeafNodes;
3330 inputTree.getNodes(inputLeafNodes);
3333 inputTree, inputLeafNodes, intersectionTree,
isovalue);
3335 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, inputLeafNodes.size()), op);
3344 template<
typename InputTreeType>
3350 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
3354 const InputTreeType& inputTree,
3355 const std::vector<BoolLeafNodeType*>& nodes,
3360 void operator()(
const tbb::blocked_range<size_t>&);
3362 mIntersectionAccessor.
tree().merge(rhs.mIntersectionAccessor.
tree());
3376 template<
typename InputTreeType>
3378 const InputTreeType& inputTree,
3379 const std::vector<BoolLeafNodeType*>& nodes,
3382 : mInputAccessor(inputTree)
3383 , mNodes(nodes.empty() ? nullptr : &nodes.front())
3384 , mIntersectionTree(false)
3385 , mIntersectionAccessor(intersectionTree)
3391 template<
typename InputTreeType>
3394 : mInputAccessor(rhs.mInputAccessor.tree())
3395 , mNodes(rhs.mNodes)
3396 , mIntersectionTree(false)
3397 , mIntersectionAccessor(mIntersectionTree)
3398 , mIsovalue(rhs.mIsovalue)
3403 template<
typename InputTreeType>
3414 for (
size_t n = range.begin(); n != range.end(); ++n) {
3418 for (
typename BoolLeafNodeType::ValueOnCIter it = node.cbeginValueOn(); it; ++it) {
3420 if (!it.getValue()) {
3422 ijk = it.getCoord();
3424 const bool inside =
isInsideValue(mInputAccessor.getValue(ijk), iso);
3443 template<
typename BoolTreeType>
3449 const std::vector<BoolLeafNodeType*>& maskNodes,
3450 BoolTreeType& borderTree)
3451 : mMaskTree(&maskTree)
3452 , mMaskNodes(maskNodes.empty() ? nullptr : &maskNodes.front())
3453 , mTmpBorderTree(false)
3454 , mBorderTree(&borderTree)
3459 : mMaskTree(rhs.mMaskTree)
3460 , mMaskNodes(rhs.mMaskNodes)
3461 , mTmpBorderTree(false)
3462 , mBorderTree(&mTmpBorderTree)
3474 for (
size_t n = range.begin(); n != range.end(); ++n) {
3478 for (
typename BoolLeafNodeType::ValueOnCIter it = node.cbeginValueOn(); it; ++it) {
3480 ijk = it.getCoord();
3482 const bool lhs = it.getValue();
3485 bool isEdgeVoxel =
false;
3488 isEdgeVoxel = (maskAcc.
probeValue(ijk, rhs) && lhs != rhs);
3491 isEdgeVoxel = isEdgeVoxel || (maskAcc.
probeValue(ijk, rhs) && lhs != rhs);
3494 isEdgeVoxel = isEdgeVoxel || (maskAcc.
probeValue(ijk, rhs) && lhs != rhs);
3497 isEdgeVoxel = isEdgeVoxel || (maskAcc.
probeValue(ijk, rhs) && lhs != rhs);
3501 isEdgeVoxel = isEdgeVoxel || (maskAcc.
probeValue(ijk, rhs) && lhs != rhs);
3504 isEdgeVoxel = isEdgeVoxel || (maskAcc.
probeValue(ijk, rhs) && lhs != rhs);
3507 isEdgeVoxel = isEdgeVoxel || (maskAcc.
probeValue(ijk, rhs) && lhs != rhs);
3518 BoolTreeType
const *
const mMaskTree;
3521 BoolTreeType mTmpBorderTree;
3522 BoolTreeType *
const mBorderTree;
3526 template<
typename BoolTreeType>
3531 SyncMaskValues(
const std::vector<BoolLeafNodeType*>& nodes,
const BoolTreeType& mask)
3532 : mNodes(nodes.empty() ? nullptr : &nodes.front())
3539 using ValueOnIter =
typename BoolLeafNodeType::ValueOnIter;
3543 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
3550 for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3551 const Index pos = it.pos();
3552 if (maskNode->getValue(pos)) {
3553 node.setValueOnly(pos,
true);
3562 BoolTreeType
const *
const mMaskTree;
3569 template<
typename BoolTreeType>
3574 MaskSurface(
const std::vector<BoolLeafNodeType*>& nodes,
const BoolTreeType& mask,
3576 : mNodes(nodes.empty() ? nullptr : &nodes.front())
3578 , mInputTransform(inputTransform)
3579 , mMaskTransform(maskTransform)
3580 , mInvertMask(invert)
3586 using ValueOnIter =
typename BoolLeafNodeType::ValueOnIter;
3590 const bool matchingTransforms = mInputTransform == mMaskTransform;
3591 const bool maskState = mInvertMask;
3593 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
3597 if (matchingTransforms) {
3603 for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3604 const Index pos = it.pos();
3605 if (maskNode->isValueOn(pos) == maskState) {
3606 node.setValueOnly(pos,
true);
3612 if (maskTreeAcc.
isValueOn(node.origin()) == maskState) {
3613 for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3614 node.setValueOnly(it.pos(),
true);
3624 for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3629 if (maskTreeAcc.
isValueOn(ijk) == maskState) {
3630 node.setValueOnly(it.pos(),
true);
3640 BoolTreeType
const *
const mMaskTree;
3643 bool const mInvertMask;
3647 template<
typename InputGr
idType>
3650 typename InputGridType::TreeType::template ValueConverter<bool>::Type& intersectionTree,
3651 typename InputGridType::TreeType::template ValueConverter<bool>::Type& borderTree,
3652 const InputGridType& inputGrid,
3655 typename InputGridType::ValueType
isovalue)
3657 using InputTreeType =
typename InputGridType::TreeType;
3658 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
3659 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
3662 if (maskGrid && maskGrid->type() == BoolGridType::gridType()) {
3665 const InputTreeType& inputTree = inputGrid.tree();
3667 const BoolGridType * surfaceMask =
static_cast<const BoolGridType*
>(maskGrid.get());
3669 const BoolTreeType& maskTree = surfaceMask->tree();
3675 std::vector<BoolLeafNodeType*> intersectionLeafNodes;
3676 intersectionTree.getNodes(intersectionLeafNodes);
3678 tbb::parallel_for(tbb::blocked_range<size_t>(0, intersectionLeafNodes.size()),
3680 intersectionLeafNodes, maskTree, transform, maskTransform, invertMask));
3686 intersectionTree, intersectionLeafNodes, borderTree);
3687 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, intersectionLeafNodes.size()), borderOp);
3692 BoolTreeType tmpIntersectionTree(
false);
3695 inputTree, intersectionLeafNodes, tmpIntersectionTree,
isovalue);
3697 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, intersectionLeafNodes.size()), op);
3699 std::vector<BoolLeafNodeType*> tmpIntersectionLeafNodes;
3700 tmpIntersectionTree.getNodes(tmpIntersectionLeafNodes);
3702 tbb::parallel_for(tbb::blocked_range<size_t>(0, tmpIntersectionLeafNodes.size()),
3705 intersectionTree.clear();
3706 intersectionTree.merge(tmpIntersectionTree);
3714 template<
typename InputTreeType>
3722 using Int16TreeType =
typename InputTreeType::template ValueConverter<Int16>::Type;
3727 const std::vector<const BoolLeafNodeType*>& intersectionLeafNodes,
3733 void operator()(
const tbb::blocked_range<size_t>&);
3735 mSignFlagsAccessor.
tree().merge(rhs.mSignFlagsAccessor.
tree());
3736 mPointIndexAccessor.
tree().merge(rhs.mPointIndexAccessor.
tree());
3752 template<
typename InputTreeType>
3754 const InputTreeType& inputTree,
3755 const std::vector<const BoolLeafNodeType*>& intersectionLeafNodes,
3759 : mInputAccessor(inputTree)
3760 , mIntersectionNodes(&intersectionLeafNodes.front())
3762 , mSignFlagsAccessor(signFlagsTree)
3764 , mPointIndexAccessor(pointIndexTree)
3771 template<
typename InputTreeType>
3773 : mInputAccessor(rhs.mInputAccessor.tree())
3774 , mIntersectionNodes(rhs.mIntersectionNodes)
3776 , mSignFlagsAccessor(mSignFlagsTree)
3778 , mPointIndexAccessor(mPointIndexTree)
3779 , mIsovalue(rhs.mIsovalue)
3784 template<
typename InputTreeType>
3788 using Int16LeafNodeType =
typename Int16TreeType::LeafNodeType;
3792 typename std::unique_ptr<Int16LeafNodeType> signsNodePt(
new Int16LeafNodeType(ijk, 0));
3794 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
3801 if (!signsNodePt.get()) signsNodePt.
reset(
new Int16LeafNodeType(origin, 0));
3802 else signsNodePt->setOrigin(origin);
3804 bool updatedNode =
false;
3808 const Index pos = it.pos();
3809 ijk = BoolLeafNodeType::offsetToLocalCoord(pos);
3812 ijk[0] <
int(BoolLeafNodeType::DIM - 1) &&
3813 ijk[1] <
int(BoolLeafNodeType::DIM - 1) &&
3814 ijk[2] <
int(BoolLeafNodeType::DIM - 1) ) {
3822 if (signFlags != 0 && signFlags != 0xFF) {
3824 const bool inside = signFlags & 0x1;
3826 int edgeFlags = inside ?
INSIDE : 0;
3828 if (!it.getValue()) {
3829 edgeFlags |= inside != ((signFlags & 0x02) != 0) ?
XEDGE : 0;
3830 edgeFlags |= inside != ((signFlags & 0x10) != 0) ?
YEDGE : 0;
3831 edgeFlags |= inside != ((signFlags & 0x08) != 0) ?
ZEDGE : 0;
3835 if (ambiguousCellFlags != 0) {
3837 origin + ijk, mIsovalue);
3840 edgeFlags |= int(signFlags);
3842 signsNodePt->setValueOn(pos,
Int16(edgeFlags));
3848 typename Index32TreeType::LeafNodeType* idxNode = mPointIndexAccessor.touchLeaf(origin);
3849 idxNode->topologyUnion(*signsNodePt);
3852 for (
auto it = idxNode->beginValueOn(); it; ++it) {
3853 idxNode->setValueOnly(it.pos(), 0);
3856 mSignFlagsAccessor.addLeaf(signsNodePt.release());
3862 template<
typename InputTreeType>
3865 typename InputTreeType::template ValueConverter<Int16>::Type& signFlagsTree,
3866 typename InputTreeType::template ValueConverter<Index32>::Type& pointIndexTree,
3867 const typename InputTreeType::template ValueConverter<bool>::Type& intersectionTree,
3868 const InputTreeType& inputTree,
3869 typename InputTreeType::ValueType
isovalue)
3871 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
3872 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
3874 std::vector<const BoolLeafNodeType*> intersectionLeafNodes;
3875 intersectionTree.getNodes(intersectionLeafNodes);
3878 inputTree, intersectionLeafNodes, signFlagsTree, pointIndexTree,
isovalue);
3880 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, intersectionLeafNodes.size()), op);
3887 template<Index32 LeafNodeLog2Dim>
3893 std::unique_ptr<
Index32[]>& leafNodeCount)
3894 : mLeafNodes(leafNodes.empty() ? nullptr : &leafNodes.front())
3895 , mData(leafNodeCount.get())
3901 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
3923 template<
typename Po
intIndexLeafNode>
3929 const std::vector<Int16LeafNodeType*>& signDataNodes,
3930 std::unique_ptr<
Index32[]>& leafNodeCount)
3931 : mPointIndexNodes(pointIndexNodes.empty() ? nullptr : &pointIndexNodes.front())
3932 , mSignDataNodes(signDataNodes.empty() ? nullptr : &signDataNodes.front())
3933 , mData(leafNodeCount.get())
3941 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
3948 std::set<IndexType> uniqueRegions;
3957 uniqueRegions.insert(
id);
3961 mData[n] =
Index32(count + uniqueRegions.size());
3972 template<
typename Po
intIndexLeafNode>
3977 MapPoints(
const std::vector<PointIndexLeafNode*>& pointIndexNodes,
3978 const std::vector<Int16LeafNodeType*>& signDataNodes,
3979 std::unique_ptr<
Index32[]>& leafNodeCount)
3980 : mPointIndexNodes(pointIndexNodes.empty() ? nullptr : &pointIndexNodes.front())
3981 , mSignDataNodes(signDataNodes.empty() ? nullptr : &signDataNodes.front())
3982 , mData(leafNodeCount.get())
3988 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
3993 Index32 pointOffset = mData[n];
3996 const Index pos = it.pos();
4013 template<
typename TreeType,
typename PrimBuilder>
4024 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
4028 bool invertSurfaceOrientation);
4032 void operator()(
const tbb::blocked_range<size_t>&)
const;
4040 bool const mInvertSurfaceOrientation;
4044 template<
typename TreeType,
typename PrimBuilder>
4046 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
4050 bool invertSurfaceOrientation)
4051 : mSignFlagsLeafNodes(signFlagsLeafNodes.empty() ? nullptr : &signFlagsLeafNodes.front())
4052 , mSignFlagsTree(&signFlagsTree)
4053 , mRefSignFlagsTree(nullptr)
4054 , mIndexTree(&idxTree)
4055 , mPolygonPoolList(&polygons)
4056 , mInvertSurfaceOrientation(invertSurfaceOrientation)
4060 template<
typename InputTreeType,
typename PrimBuilder>
4065 Int16ValueAccessor signAcc(*mSignFlagsTree);
4069 const bool invertSurfaceOrientation = mInvertSurfaceOrientation;
4076 std::unique_ptr<Int16ValueAccessor> refSignAcc;
4077 if (mRefSignFlagsTree) refSignAcc.
reset(
new Int16ValueAccessor(*mRefSignFlagsTree));
4079 for (
size_t n = range.begin(); n != range.end(); ++n) {
4082 origin = node.origin();
4086 typename Int16LeafNodeType::ValueOnCIter iter = node.cbeginValueOn();
4087 for (; iter; ++iter) {
4088 if (iter.getValue() &
XEDGE) ++edgeCount;
4089 if (iter.getValue() &
YEDGE) ++edgeCount;
4090 if (iter.getValue() &
ZEDGE) ++edgeCount;
4093 if(edgeCount == 0)
continue;
4095 mesher.init(edgeCount, (*mPolygonPoolList)[n]);
4100 if (!signleafPt || !idxLeafPt)
continue;
4104 if (refSignAcc) refSignLeafPt = refSignAcc->probeConstLeaf(origin);
4108 for (iter = node.cbeginValueOn(); iter; ++iter) {
4109 ijk = iter.getCoord();
4111 Int16 flags = iter.getValue();
4113 if (!(flags & 0xE00))
continue;
4116 if (refSignLeafPt) {
4117 refFlags = refSignLeafPt->getValue(iter.pos());
4124 const uint8_t cell = uint8_t(
SIGNS & flags);
4132 if (ijk[0] > origin[0] && ijk[1] > origin[1] && ijk[2] > origin[2]) {
4134 flags, refFlags, offsets, ijk, *signleafPt, *idxLeafPt, mesher);
4137 flags, refFlags, offsets, ijk, signAcc, idxAcc, mesher);
4150 template<
typename T>
4153 CopyArray(T * outputArray,
const T * inputArray,
size_t outputOffset = 0)
4154 : mOutputArray(outputArray), mInputArray(inputArray), mOutputOffset(outputOffset)
4158 void operator()(
const tbb::blocked_range<size_t>& inputArrayRange)
const
4160 const size_t offset = mOutputOffset;
4161 for (
size_t n = inputArrayRange.begin(), N = inputArrayRange.end(); n < N; ++n) {
4162 mOutputArray[offset + n] = mInputArray[n];
4167 T *
const mOutputArray;
4168 T
const *
const mInputArray;
4169 size_t const mOutputOffset;
4176 const std::vector<uint8_t>& pointFlags,
4178 std::unique_ptr<
unsigned[]>& numQuadsToDivide)
4179 : mPolygonPoolList(&polygons)
4180 , mPointFlags(pointFlags.empty() ? nullptr : &pointFlags.front())
4181 , mPoints(points.get())
4182 , mNumQuadsToDivide(numQuadsToDivide.get())
4188 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
4195 for (
size_t i = 0, I = polygons.
numQuads(); i < I; ++i) {
4203 const bool edgePoly = mPointFlags[quad[0]] || mPointFlags[quad[1]]
4204 || mPointFlags[quad[2]] || mPointFlags[quad[3]];
4206 if (!edgePoly)
continue;
4208 const Vec3s& p0 = mPoints[quad[0]];
4209 const Vec3s& p1 = mPoints[quad[1]];
4210 const Vec3s& p2 = mPoints[quad[2]];
4211 const Vec3s& p3 = mPoints[quad[3]];
4220 mNumQuadsToDivide[n] = count;
4226 uint8_t
const *
const mPointFlags;
4227 Vec3s const *
const mPoints;
4228 unsigned *
const mNumQuadsToDivide;
4238 std::unique_ptr<
unsigned[]>& numQuadsToDivide,
4239 std::unique_ptr<
unsigned[]>& centroidOffsets)
4240 : mPolygonPoolList(&polygons)
4241 , mPoints(points.get())
4242 , mCentroids(centroids.get())
4243 , mNumQuadsToDivide(numQuadsToDivide.get())
4244 , mCentroidOffsets(centroidOffsets.get())
4251 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
4255 const size_t nonplanarCount = size_t(mNumQuadsToDivide[n]);
4257 if (nonplanarCount > 0) {
4263 size_t offset = mCentroidOffsets[n];
4265 size_t triangleIdx = 0;
4267 for (
size_t i = 0, I = polygons.
numQuads(); i < I; ++i) {
4269 const char quadFlags = polygons.
quadFlags(i);
4272 unsigned newPointIdx = unsigned(offset + mPointCount);
4276 mCentroids[offset] = (mPoints[quad[0]] + mPoints[quad[1]] +
4277 mPoints[quad[2]] + mPoints[quad[3]]) * 0.25f;
4284 triangle[0] = quad[0];
4285 triangle[1] = newPointIdx;
4286 triangle[2] = quad[3];
4296 triangle[0] = quad[0];
4297 triangle[1] = quad[1];
4298 triangle[2] = newPointIdx;
4308 triangle[0] = quad[1];
4309 triangle[1] = quad[2];
4310 triangle[2] = newPointIdx;
4321 triangle[0] = quad[2];
4322 triangle[1] = quad[3];
4323 triangle[2] = newPointIdx;
4334 for (
size_t i = 0, I = polygons.
numTriangles(); i < I; ++i) {
4341 for (
size_t i = 0, I = polygons.
numQuads(); i < I; ++i) {
4345 tmpPolygons.
quad(quadIdx) = quad;
4351 polygons.
copy(tmpPolygons);
4358 Vec3s const *
const mPoints;
4359 Vec3s *
const mCentroids;
4360 unsigned *
const mNumQuadsToDivide;
4361 unsigned *
const mCentroidOffsets;
4362 size_t const mPointCount;
4369 size_t polygonPoolListSize,
4371 size_t& pointListSize,
4372 std::vector<uint8_t>& pointFlags)
4374 const tbb::blocked_range<size_t> polygonPoolListRange(0, polygonPoolListSize);
4376 std::unique_ptr<unsigned[]> numQuadsToDivide(
new unsigned[polygonPoolListSize]);
4378 tbb::parallel_for(polygonPoolListRange,
4381 std::unique_ptr<unsigned[]> centroidOffsets(
new unsigned[polygonPoolListSize]);
4383 size_t centroidCount = 0;
4387 for (
size_t n = 0, N = polygonPoolListSize; n < N; ++n) {
4388 centroidOffsets[n] = sum;
4389 sum += numQuadsToDivide[n];
4391 centroidCount = size_t(sum);
4394 std::unique_ptr<Vec3s[]> centroidList(
new Vec3s[centroidCount]);
4396 tbb::parallel_for(polygonPoolListRange,
4398 centroidList, numQuadsToDivide, centroidOffsets));
4400 if (centroidCount > 0) {
4402 const size_t newPointListSize = centroidCount + pointListSize;
4404 std::unique_ptr<openvdb::Vec3s[]> newPointList(
new openvdb::Vec3s[newPointListSize]);
4406 tbb::parallel_for(tbb::blocked_range<size_t>(0, pointListSize),
4409 tbb::parallel_for(tbb::blocked_range<size_t>(0, newPointListSize - pointListSize),
4412 pointListSize = newPointListSize;
4413 pointList.swap(newPointList);
4414 pointFlags.resize(pointListSize, 0);
4422 const std::vector<uint8_t>& pointFlags)
4423 : mPolygonPoolList(&polygons)
4424 , mPointFlags(pointFlags.empty() ? nullptr : &pointFlags.front())
4430 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
4434 for (
size_t i = 0, I = polygons.
numQuads(); i < I; ++i) {
4442 const bool hasSeamLinePoint =
4443 mPointFlags[verts[0]] || mPointFlags[verts[1]] ||
4444 mPointFlags[verts[2]] || mPointFlags[verts[3]];
4446 if (!hasSeamLinePoint) {
4452 for (
size_t i = 0, I = polygons.
numTriangles(); i < I; ++i) {
4460 const bool hasSeamLinePoint =
4461 mPointFlags[verts[0]] || mPointFlags[verts[1]] || mPointFlags[verts[2]];
4463 if (!hasSeamLinePoint) {
4475 uint8_t
const *
const mPointFlags;
4481 std::vector<uint8_t>& pointFlags)
4483 tbb::parallel_for(tbb::blocked_range<size_t>(0, polygonPoolListSize),
4491 template<
typename InputTreeType>
4495 const PointList& pointList, std::unique_ptr<uint8_t[]>& pointMask,
4497 : mInputTree(&inputTree)
4498 , mPolygonPoolList(&polygons)
4499 , mPointList(&pointList)
4500 , mPointMask(pointMask.get())
4501 , mTransform(transform)
4502 , mInvertSurfaceOrientation(invertSurfaceOrientation)
4508 using ValueType =
typename InputTreeType::LeafNodeType::ValueType;
4511 Vec3s centroid, normal;
4514 const bool invertGradientDir = mInvertSurfaceOrientation || isBoolValue<ValueType>();
4516 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
4518 const PolygonPool& polygons = (*mPolygonPoolList)[n];
4520 for (
size_t i = 0, I = polygons.
numTriangles(); i < I; ++i) {
4524 const Vec3s& v0 = (*mPointList)[verts[0]];
4525 const Vec3s& v1 = (*mPointList)[verts[1]];
4526 const Vec3s& v2 = (*mPointList)[verts[2]];
4528 normal = (v2 - v0).cross((v1 - v0));
4531 centroid = (v0 + v1 + v2) * (1.0f / 3.0f);
4537 if (invertGradientDir) {
4542 if (dir.dot(normal) < -0.5f) {
4547 mPointMask[verts[0]] = 1;
4548 mPointMask[verts[1]] = 1;
4549 mPointMask[verts[2]] = 1;
4558 InputTreeType
const *
const mInputTree;
4561 uint8_t *
const mPointMask;
4563 bool const mInvertSurfaceOrientation;
4567 template<
typename InputTree>
4570 bool invertSurfaceOrientation,
4571 const InputTree& inputTree,
4574 size_t polygonPoolListSize,
4576 const size_t pointListSize)
4578 const tbb::blocked_range<size_t> polygonPoolListRange(0, polygonPoolListSize);
4580 std::unique_ptr<uint8_t[]> pointMask(
new uint8_t[pointListSize]);
4581 fillArray(pointMask.get(), uint8_t(0), pointListSize);
4583 tbb::parallel_for(polygonPoolListRange,
4585 inputTree, polygonPoolList, pointList, pointMask, transform, invertSurfaceOrientation));
4587 std::unique_ptr<uint8_t[]> pointUpdates(
new uint8_t[pointListSize]);
4588 fillArray(pointUpdates.get(), uint8_t(0), pointListSize);
4590 std::unique_ptr<Vec3s[]> newPoints(
new Vec3s[pointListSize]);
4591 fillArray(newPoints.get(),
Vec3s(0.0f, 0.0f, 0.0f), pointListSize);
4593 for (
size_t n = 0, N = polygonPoolListSize; n < N; ++n) {
4597 for (
size_t i = 0, I = polygons.
numQuads(); i < I; ++i) {
4600 for (
int v = 0; v < 4; ++v) {
4602 const unsigned pointIdx = verts[v];
4604 if (pointMask[pointIdx] == 1) {
4606 newPoints[pointIdx] +=
4607 pointList[verts[0]] + pointList[verts[1]] +
4608 pointList[verts[2]] + pointList[verts[3]];
4610 pointUpdates[pointIdx] = uint8_t(pointUpdates[pointIdx] + 4);
4615 for (
size_t i = 0, I = polygons.
numTriangles(); i < I; ++i) {
4618 for (
int v = 0; v < 3; ++v) {
4620 const unsigned pointIdx = verts[v];
4622 if (pointMask[pointIdx] == 1) {
4623 newPoints[pointIdx] +=
4624 pointList[verts[0]] + pointList[verts[1]] + pointList[verts[2]];
4626 pointUpdates[pointIdx] = uint8_t(pointUpdates[pointIdx] + 3);
4632 for (
size_t n = 0, N = pointListSize; n < N; ++n) {
4633 if (pointUpdates[n] > 0) {
4634 const double weight = 1.0 / double(pointUpdates[n]);
4635 pointList[n] = newPoints[n] * float(weight);
4652 , mTriangles(nullptr)
4653 , mQuadFlags(nullptr)
4654 , mTriangleFlags(nullptr)
4661 : mNumQuads(numQuads)
4662 , mNumTriangles(numTriangles)
4665 , mQuadFlags(new char[mNumQuads])
4666 , mTriangleFlags(new char[mNumTriangles])
4677 for (
size_t i = 0; i < mNumQuads; ++i) {
4678 mQuads[i] = rhs.mQuads[i];
4679 mQuadFlags[i] = rhs.mQuadFlags[i];
4682 for (
size_t i = 0; i < mNumTriangles; ++i) {
4683 mTriangles[i] = rhs.mTriangles[i];
4684 mTriangleFlags[i] = rhs.mTriangleFlags[i];
4694 mQuadFlags.reset(
new char[mNumQuads]);
4702 mQuads.reset(
nullptr);
4703 mQuadFlags.reset(
nullptr);
4710 mNumTriangles = size;
4712 mTriangleFlags.reset(
new char[mNumTriangles]);
4720 mTriangles.reset(
nullptr);
4721 mTriangleFlags.reset(
nullptr);
4728 if (!(n < mNumQuads))
return false;
4733 mQuads.reset(
nullptr);
4737 std::unique_ptr<char[]> flags(
new char[n]);
4739 for (
size_t i = 0; i < n; ++i) {
4740 quads[i] = mQuads[i];
4741 flags[i] = mQuadFlags[i];
4745 mQuadFlags.swap(flags);
4757 if (!(n < mNumTriangles))
return false;
4762 mTriangles.reset(
nullptr);
4765 std::unique_ptr<openvdb::Vec3I[]> triangles(
new openvdb::Vec3I[n]);
4766 std::unique_ptr<char[]> flags(
new char[n]);
4768 for (
size_t i = 0; i < n; ++i) {
4769 triangles[i] = mTriangles[i];
4770 flags[i] = mTriangleFlags[i];
4773 mTriangles.swap(triangles);
4774 mTriangleFlags.swap(flags);
4791 , mSeamPointListSize(0)
4792 , mPolygonPoolListSize(0)
4794 , mPrimAdaptivity(adaptivity)
4795 , mSecAdaptivity(0.0)
4797 , mSurfaceMaskGrid(
GridBase::ConstPtr())
4798 , mAdaptivityGrid(
GridBase::ConstPtr())
4799 , mAdaptivityMaskTree(
TreeBase::ConstPtr())
4802 , mInvertSurfaceMask(false)
4804 , mQuantizedSeamPoints(nullptr)
4814 mSecAdaptivity = secAdaptivity;
4819 mSeamPointListSize = 0;
4820 mQuantizedSeamPoints.reset(
nullptr);
4827 mSurfaceMaskGrid = mask;
4828 mInvertSurfaceMask = invertMask;
4835 mAdaptivityGrid = grid;
4842 mAdaptivityMaskTree = tree;
4846 template<
typename InputGr
idType>
4852 using InputTreeType =
typename InputGridType::TreeType;
4853 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
4854 using InputValueType =
typename InputLeafNodeType::ValueType;
4858 using FloatTreeType =
typename InputTreeType::template ValueConverter<float>::Type;
4860 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
4861 using Int16TreeType =
typename InputTreeType::template ValueConverter<Int16>::Type;
4862 using Int16LeafNodeType =
typename Int16TreeType::LeafNodeType;
4863 using Index32TreeType =
typename InputTreeType::template ValueConverter<Index32>::Type;
4864 using Index32LeafNodeType =
typename Index32TreeType::LeafNodeType;
4869 mPolygonPoolListSize = 0;
4871 mPointFlags.clear();
4876 const InputValueType
isovalue = InputValueType(mIsovalue);
4877 const float adaptivityThreshold = float(mPrimAdaptivity);
4878 const bool adaptive = mPrimAdaptivity > 1e-7 || mSecAdaptivity > 1e-7;
4884 const bool invertSurfaceOrientation = (!volume_to_mesh_internal::isBoolValue<InputValueType>()
4889 const InputTreeType& inputTree = inputGrid.tree();
4891 BoolTreeType intersectionTree(
false), adaptivityMask(
false);
4893 if (mAdaptivityMaskTree && mAdaptivityMaskTree->type() == BoolTreeType::treeType()) {
4894 const BoolTreeType *refAdaptivityMask=
4895 static_cast<const BoolTreeType*
>(mAdaptivityMaskTree.get());
4896 adaptivityMask.topologyUnion(*refAdaptivityMask);
4899 Int16TreeType signFlagsTree(0);
4906 intersectionTree, inputTree,
isovalue);
4909 inputGrid, mSurfaceMaskGrid, mInvertSurfaceMask,
isovalue);
4911 if (intersectionTree.empty())
return;
4914 signFlagsTree, pointIndexTree, intersectionTree, inputTree,
isovalue);
4916 intersectionTree.clear();
4918 std::vector<Index32LeafNodeType*> pointIndexLeafNodes;
4919 pointIndexTree.getNodes(pointIndexLeafNodes);
4921 std::vector<Int16LeafNodeType*> signFlagsLeafNodes;
4922 signFlagsTree.getNodes(signFlagsLeafNodes);
4924 const tbb::blocked_range<size_t> auxiliaryLeafNodeRange(0, signFlagsLeafNodes.size());
4929 Int16TreeType* refSignFlagsTree =
nullptr;
4930 Index32TreeType* refPointIndexTree =
nullptr;
4931 InputTreeType
const* refInputTree =
nullptr;
4933 if (mRefGrid && mRefGrid->type() == InputGridType::gridType()) {
4935 const InputGridType* refGrid =
static_cast<const InputGridType*
>(mRefGrid.get());
4936 refInputTree = &refGrid->tree();
4938 if (!mRefSignTree && !mRefIdxTree) {
4942 typename Int16TreeType::Ptr refSignFlagsTreePt(
new Int16TreeType(0));
4943 typename Index32TreeType::Ptr refPointIndexTreePt(
4946 BoolTreeType refIntersectionTree(
false);
4949 refIntersectionTree, *refInputTree,
isovalue);
4952 *refPointIndexTreePt, refIntersectionTree, *refInputTree,
isovalue);
4954 mRefSignTree = refSignFlagsTreePt;
4955 mRefIdxTree = refPointIndexTreePt;
4958 if (mRefSignTree && mRefIdxTree) {
4962 refSignFlagsTree =
static_cast<Int16TreeType*
>(mRefSignTree.get());
4963 refPointIndexTree =
static_cast<Index32TreeType*
>(mRefIdxTree.get());
4967 if (refSignFlagsTree && refPointIndexTree) {
4973 if (mSeamPointListSize == 0) {
4977 std::vector<Int16LeafNodeType*> refSignFlagsLeafNodes;
4978 refSignFlagsTree->getNodes(refSignFlagsLeafNodes);
4980 std::unique_ptr<Index32[]> leafNodeOffsets(
4981 new Index32[refSignFlagsLeafNodes.size()]);
4983 tbb::parallel_for(tbb::blocked_range<size_t>(0, refSignFlagsLeafNodes.size()),
4985 refSignFlagsLeafNodes, leafNodeOffsets));
4989 for (
size_t n = 0, N = refSignFlagsLeafNodes.size(); n < N; ++n) {
4990 const Index32 tmp = leafNodeOffsets[n];
4991 leafNodeOffsets[n] = count;
4994 mSeamPointListSize = size_t(count);
4997 if (mSeamPointListSize != 0) {
4999 mQuantizedSeamPoints.reset(
new uint32_t[mSeamPointListSize]);
5001 memset(mQuantizedSeamPoints.get(), 0,
sizeof(uint32_t) * mSeamPointListSize);
5003 std::vector<Index32LeafNodeType*> refPointIndexLeafNodes;
5004 refPointIndexTree->getNodes(refPointIndexLeafNodes);
5006 tbb::parallel_for(tbb::blocked_range<size_t>(0, refPointIndexLeafNodes.size()),
5008 refPointIndexLeafNodes, refSignFlagsLeafNodes, leafNodeOffsets));
5012 if (mSeamPointListSize != 0) {
5014 tbb::parallel_for(auxiliaryLeafNodeRange,
5016 signFlagsLeafNodes, inputTree, *refPointIndexTree, *refSignFlagsTree,
5017 mQuantizedSeamPoints.get(),
isovalue));
5022 const bool referenceMeshing = refSignFlagsTree && refPointIndexTree && refInputTree;
5027 std::unique_ptr<Index32[]> leafNodeOffsets(
new Index32[signFlagsLeafNodes.size()]);
5031 inputGrid, pointIndexTree, pointIndexLeafNodes, signFlagsLeafNodes,
5032 isovalue, adaptivityThreshold, invertSurfaceOrientation);
5034 if (mAdaptivityGrid && mAdaptivityGrid->type() == FloatGridType::gridType()) {
5035 const FloatGridType* adaptivityGrid =
5036 static_cast<const FloatGridType*
>(mAdaptivityGrid.get());
5040 if (!adaptivityMask.empty()) {
5044 if (referenceMeshing) {
5048 tbb::parallel_for(auxiliaryLeafNodeRange, mergeOp);
5051 op(pointIndexLeafNodes, signFlagsLeafNodes, leafNodeOffsets);
5053 tbb::parallel_for(auxiliaryLeafNodeRange, op);
5058 op(signFlagsLeafNodes, leafNodeOffsets);
5060 tbb::parallel_for(auxiliaryLeafNodeRange, op);
5066 for (
size_t n = 0, N = signFlagsLeafNodes.size(); n < N; ++n) {
5067 const Index32 tmp = leafNodeOffsets[n];
5074 mPointFlags.clear();
5082 op(mPoints.get(), inputTree, pointIndexLeafNodes,
5083 signFlagsLeafNodes, leafNodeOffsets, transform, mIsovalue);
5085 if (referenceMeshing) {
5086 mPointFlags.resize(mPointListSize);
5087 op.
setRefData(*refInputTree, *refPointIndexTree, *refSignFlagsTree,
5088 mQuantizedSeamPoints.get(), &mPointFlags.front());
5091 tbb::parallel_for(auxiliaryLeafNodeRange, op);
5097 mPolygonPoolListSize = signFlagsLeafNodes.size();
5098 mPolygons.reset(
new PolygonPool[mPolygonPoolListSize]);
5105 op(signFlagsLeafNodes, signFlagsTree, pointIndexTree,
5106 mPolygons, invertSurfaceOrientation);
5108 if (referenceMeshing) {
5112 tbb::parallel_for(auxiliaryLeafNodeRange, op);
5119 op(signFlagsLeafNodes, signFlagsTree, pointIndexTree,
5120 mPolygons, invertSurfaceOrientation);
5122 if (referenceMeshing) {
5126 tbb::parallel_for(auxiliaryLeafNodeRange, op);
5130 signFlagsTree.clear();
5131 pointIndexTree.clear();
5134 if (adaptive && mRelaxDisorientedTriangles) {
5136 inputTree, transform, mPolygons, mPolygonPoolListSize, mPoints, mPointListSize);
5140 if (referenceMeshing) {
5142 mPolygons, mPolygonPoolListSize, mPoints, mPointListSize, mPointFlags);
5157 template<
typename Gr
idType>
5158 inline typename std::enable_if<std::is_scalar<typename GridType::ValueType>::value,
void>
::type
5160 const GridType& grid,
5161 std::vector<Vec3s>& points,
5162 std::vector<Vec3I>& triangles,
5163 std::vector<Vec4I>& quads,
5173 points.resize(mesher.pointListSize());
5177 tbb::parallel_for(tbb::blocked_range<size_t>(0, points.size()), ptnCpy);
5178 mesher.pointList().reset(
nullptr);
5184 size_t numQuads = 0, numTriangles = 0;
5185 for (
size_t n = 0, N = mesher.polygonPoolListSize(); n < N; ++n) {
5186 openvdb::tools::PolygonPool& polygons = polygonPoolList[n];
5187 numTriangles += polygons.numTriangles();
5188 numQuads += polygons.numQuads();
5192 triangles.resize(numTriangles);
5194 quads.resize(numQuads);
5198 size_t qIdx = 0, tIdx = 0;
5199 for (
size_t n = 0, N = mesher.polygonPoolListSize(); n < N; ++n) {
5200 openvdb::tools::PolygonPool& polygons = polygonPoolList[n];
5202 for (
size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
5203 quads[qIdx++] = polygons.quad(i);
5206 for (
size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
5207 triangles[tIdx++] = polygons.triangle(i);
5213 template<
typename Gr
idType>
5214 inline typename std::enable_if<!std::is_scalar<typename GridType::ValueType>::value,
void>
::type
5217 std::vector<Vec3s>&,
5218 std::vector<Vec3I>&,
5219 std::vector<Vec4I>&,
5224 OPENVDB_THROW(TypeError,
"volume to mesh conversion is supported only for scalar grids");
5231 template<
typename Gr
idType>
5234 const GridType& grid,
5235 std::vector<Vec3s>& points,
5236 std::vector<Vec3I>& triangles,
5237 std::vector<Vec4I>& quads,
5246 template<
typename Gr
idType>
5249 const GridType& grid,
5250 std::vector<Vec3s>& points,
5251 std::vector<Vec4I>& quads,
5254 std::vector<Vec3I> triangles;
5255 doVolumeToMesh(grid, points, triangles, quads,
isovalue, 0.0,
true);
5262 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED