10 #ifndef OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED 11 #define OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED 20 #include <tbb/blocked_range.h> 21 #include <tbb/parallel_for.h> 22 #include <tbb/parallel_reduce.h> 23 #include <tbb/task_arena.h> 30 #include <type_traits> 54 template<
typename Gr
idType>
58 std::vector<Vec3s>& points,
59 std::vector<Vec4I>& quads,
60 double isovalue = 0.0);
83 template<
typename Gr
idType>
87 std::vector<Vec3s>& points,
88 std::vector<Vec3I>& triangles,
89 std::vector<Vec4I>& quads,
90 double isovalue = 0.0,
91 double adaptivity = 0.0,
92 bool relaxDisorientedTriangles =
true);
108 inline PolygonPool(
const size_t numQuads,
const size_t numTriangles);
112 inline void resetQuads(
size_t size);
113 inline void clearQuads();
115 inline void resetTriangles(
size_t size);
116 inline void clearTriangles();
121 const size_t&
numQuads()
const {
return mNumQuads; }
136 const char&
quadFlags(
size_t n)
const {
return mQuadFlags[n]; }
145 inline bool trimQuads(
const size_t n,
bool reallocate =
false);
146 inline bool trimTrinagles(
const size_t n,
bool reallocate =
false);
152 size_t mNumQuads, mNumTriangles;
153 std::unique_ptr<openvdb::Vec4I[]> mQuads;
154 std::unique_ptr<openvdb::Vec3I[]> mTriangles;
155 std::unique_ptr<char[]> mQuadFlags, mTriangleFlags;
177 VolumeToMesh(
double isovalue = 0,
double adaptivity = 0,
bool relaxDisorientedTriangles =
true);
193 const std::vector<uint8_t>&
pointFlags()
const {
return mPointFlags; }
202 template<
typename InputGr
idType>
203 void operator()(
const InputGridType&);
256 size_t mPointListSize, mSeamPointListSize, mPolygonPoolListSize;
257 double mIsovalue, mPrimAdaptivity, mSecAdaptivity;
264 bool mInvertSurfaceMask, mRelaxDisorientedTriangles;
266 std::unique_ptr<uint32_t[]> mQuantizedSeamPoints;
267 std::vector<uint8_t> mPointFlags;
281 const std::vector<Vec3d>& points,
282 const std::vector<Vec3d>& normals)
288 if (points.empty())
return avgPos;
290 for (
size_t n = 0, N = points.size(); n < N; ++n) {
294 avgPos /= double(points.size());
298 double m00=0,m01=0,m02=0,
305 for (
size_t n = 0, N = points.size(); n < N; ++n) {
307 const Vec3d& n_ref = normals[n];
310 m00 += n_ref[0] * n_ref[0];
311 m11 += n_ref[1] * n_ref[1];
312 m22 += n_ref[2] * n_ref[2];
314 m01 += n_ref[0] * n_ref[1];
315 m02 += n_ref[0] * n_ref[2];
316 m12 += n_ref[1] * n_ref[2];
319 rhs += n_ref * n_ref.
dot(points[n] - avgPos);
344 Mat3d D = Mat3d::identity();
347 double tolerance =
std::max(std::abs(eigenValues[0]), std::abs(eigenValues[1]));
348 tolerance =
std::max(tolerance, std::abs(eigenValues[2]));
352 for (
int i = 0; i < 3; ++i ) {
353 if (std::abs(eigenValues[i]) < tolerance) {
357 D[i][i] = 1.0 / eigenValues[i];
364 return avgPos + pseudoInv * rhs;
379 namespace volume_to_mesh_internal {
381 template<
typename ValueType>
384 FillArray(ValueType* array,
const ValueType& v) : mArray(array), mValue(v) { }
386 void operator()(
const tbb::blocked_range<size_t>& range)
const {
387 const ValueType v = mValue;
388 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
393 ValueType *
const mArray;
394 const ValueType mValue;
398 template<
typename ValueType>
400 fillArray(ValueType* array,
const ValueType& val,
const size_t length)
402 const auto grainSize = std::max<size_t>(
403 length / tbb::this_task_arena::max_concurrency(), 1024);
404 const tbb::blocked_range<size_t> range(0, length, grainSize);
405 tbb::parallel_for(range, FillArray<ValueType>(array, val), tbb::simple_partitioner());
410 enum { SIGNS = 0xFF, EDGES = 0xE00, INSIDE = 0x100,
411 XEDGE = 0x200, YEDGE = 0x400, ZEDGE = 0x800, SEAM = 0x1000};
415 const bool sAdaptable[256] = {
416 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,
417 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,
418 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,
419 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,
420 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,
421 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,
422 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,
423 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};
427 const unsigned char sAmbiguousFace[256] = {
428 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,
429 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,
430 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,
431 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,
432 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,
433 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,
434 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,
435 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};
441 const unsigned char sEdgeGroupTable[256][13] = {
442 {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},
443 {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},
444 {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},
445 {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},
446 {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},
447 {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},
448 {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},
449 {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},
450 {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},
451 {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},
452 {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},
453 {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},
454 {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},
455 {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},
456 {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},
457 {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},
458 {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},
459 {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},
460 {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},
461 {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},
462 {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},
463 {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},
464 {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},
465 {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},
466 {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},
467 {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},
468 {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},
469 {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},
470 {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},
471 {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},
472 {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},
473 {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},
474 {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},
475 {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},
476 {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},
477 {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},
478 {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},
479 {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},
480 {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},
481 {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},
482 {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},
483 {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},
484 {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},
485 {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},
486 {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},
487 {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},
488 {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},
489 {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},
490 {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},
491 {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},
492 {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},
493 {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},
494 {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},
495 {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},
496 {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},
497 {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},
498 {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},
499 {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},
500 {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},
501 {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},
502 {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},
503 {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},
504 {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},
505 {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},
506 {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},
507 {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},
508 {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},
509 {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},
510 {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},
511 {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},
512 {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},
513 {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},
514 {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},
515 {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},
516 {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},
517 {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},
518 {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},
519 {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},
520 {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},
521 {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},
522 {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},
523 {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},
524 {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},
525 {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},
526 {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},
527 {0,0,0,0,0,0,0,0,0,0,0,0,0}};
533 const Vec3d& p0,
const Vec3d& p1,
534 const Vec3d& p2,
const Vec3d& p3,
535 const double epsilon = 0.001)
538 Vec3d normal = (p2-p0).cross(p1-p3);
540 const Vec3d centroid = (p0 + p1 + p2 + p3);
541 const double d = centroid.
dot(normal) * 0.25;
545 double absDist = std::abs(p0.dot(normal) - d);
546 if (absDist > epsilon)
return false;
548 absDist = std::abs(p1.dot(normal) - d);
549 if (absDist > epsilon)
return false;
551 absDist = std::abs(p2.dot(normal) - d);
552 if (absDist > epsilon)
return false;
554 absDist = std::abs(p3.dot(normal) - d);
555 if (absDist > epsilon)
return false;
568 MASK_FIRST_10_BITS = 0x000003FF,
569 MASK_DIRTY_BIT = 0x80000000,
570 MASK_INVALID_BIT = 0x40000000
574 packPoint(
const Vec3d& v)
579 OPENVDB_ASSERT(!(v.x() > 1.0) && !(v.y() > 1.0) && !(v.z() > 1.0));
580 OPENVDB_ASSERT(!(v.x() < 0.0) && !(v.y() < 0.0) && !(v.z() < 0.0));
582 data |= (uint32_t(v.x() * 1023.0) & MASK_FIRST_10_BITS) << 20;
583 data |= (uint32_t(v.y() * 1023.0) & MASK_FIRST_10_BITS) << 10;
584 data |= (uint32_t(v.z() * 1023.0) & MASK_FIRST_10_BITS);
590 unpackPoint(uint32_t data)
593 v.
z() = double(data & MASK_FIRST_10_BITS) * 0.0009775171;
595 v.y() = double(data & MASK_FIRST_10_BITS) * 0.0009775171;
597 v.x() = double(data & MASK_FIRST_10_BITS) * 0.0009775171;
607 inline bool isBoolValue() {
return false; }
610 inline bool isBoolValue<bool>() {
return true; }
613 inline bool isInsideValue(T value, T isovalue) {
return value < isovalue; }
616 inline bool isInsideValue<bool>(
bool value,
bool ) {
return value; }
621 template <
typename LeafT,
622 bool IsBool = std::is_same<typename LeafT::ValueType, bool>::value>
623 struct LeafBufferAccessor
625 using T =
typename LeafT::ValueType;
626 LeafBufferAccessor(
const LeafT& leaf) : mData(leaf.buffer().data()) {}
627 inline T
get(
const Index idx)
const {
return mData[idx]; }
628 const T*
const mData;
631 template <
typename LeafT>
632 struct LeafBufferAccessor<LeafT, true>
635 LeafBufferAccessor(
const LeafT& leaf) : mLeaf(leaf) {}
636 inline T
get(
const Index idx)
const {
return mLeaf.getValue(idx); }
642 template <
typename LeafT>
643 bool isInternalLeafCoord(
const Coord& ijk)
646 ijk[0] < int(LeafT::DIM - 1) &&
647 ijk[1] < int(LeafT::DIM - 1) &&
648 ijk[2] < int(LeafT::DIM - 1);
653 template<
typename AccessorT,
typename ValueT>
655 getCellVertexValues(
const AccessorT& accessor,
657 std::array<ValueT, 8>& values)
659 values[0] = ValueT(accessor.getValue(ijk));
661 values[1] = ValueT(accessor.getValue(ijk));
663 values[2] = ValueT(accessor.getValue(ijk));
665 values[3] = ValueT(accessor.getValue(ijk));
667 values[4] = ValueT(accessor.getValue(ijk));
669 values[5] = ValueT(accessor.getValue(ijk));
671 values[6] = ValueT(accessor.getValue(ijk));
673 values[7] = ValueT(accessor.getValue(ijk));
678 template<
typename LeafT,
typename ValueT>
680 getCellVertexValues(
const LeafT& leaf,
682 std::array<ValueT, 8>& values)
684 const LeafBufferAccessor<LeafT> acc(leaf);
686 values[0] = ValueT(acc.get(offset));
687 values[3] = ValueT(acc.get(offset + 1));
688 values[4] = ValueT(acc.get(offset + LeafT::DIM));
689 values[7] = ValueT(acc.get(offset + LeafT::DIM + 1));
690 values[1] = ValueT(acc.get(offset + (LeafT::DIM * LeafT::DIM)));
691 values[2] = ValueT(acc.get(offset + (LeafT::DIM * LeafT::DIM) + 1));
692 values[5] = ValueT(acc.get(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM));
693 values[6] = ValueT(acc.get(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM + 1));
697 template<
typename ValueType>
699 computeSignFlags(
const std::array<ValueType, 8>& values,
const ValueType iso)
702 signs |= isInsideValue(values[0], iso) ? 1u : 0u;
703 signs |= isInsideValue(values[1], iso) ? 2u : 0u;
704 signs |= isInsideValue(values[2], iso) ? 4u : 0u;
705 signs |= isInsideValue(values[3], iso) ? 8u : 0u;
706 signs |= isInsideValue(values[4], iso) ? 16u : 0u;
707 signs |= isInsideValue(values[5], iso) ? 32u : 0u;
708 signs |= isInsideValue(values[6], iso) ? 64u : 0u;
709 signs |= isInsideValue(values[7], iso) ? 128u : 0u;
710 return uint8_t(signs);
716 template<
typename AccessorT>
718 evalCellSigns(
const AccessorT& accessor,
const Coord& ijk,
typename AccessorT::ValueType iso)
722 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 1u;
724 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 2u;
726 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 4u;
728 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 8u;
729 coord[1] += 1; coord[2] = ijk[2];
730 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 16u;
732 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 32u;
734 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 64u;
736 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 128u;
737 return uint8_t(signs);
743 template<
typename LeafT>
745 evalCellSigns(
const LeafT& leaf,
const Index offset,
typename LeafT::ValueType iso)
747 const LeafBufferAccessor<LeafT> acc(leaf);
750 if (isInsideValue(acc.get(offset), iso)) signs |= 1u;
751 if (isInsideValue(acc.get(offset + (LeafT::DIM * LeafT::DIM)), iso)) signs |= 2u;
752 if (isInsideValue(acc.get(offset + (LeafT::DIM * LeafT::DIM) + 1), iso)) signs |= 4u;
753 if (isInsideValue(acc.get(offset + 1), iso)) signs |= 8u;
754 if (isInsideValue(acc.get(offset + LeafT::DIM), iso)) signs |= 16u;
755 if (isInsideValue(acc.get(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM), iso)) signs |= 32u;
756 if (isInsideValue(acc.get(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM + 1), iso)) signs |= 64u;
757 if (isInsideValue(acc.get(offset + LeafT::DIM + 1), iso)) signs |= 128u;
758 return uint8_t(signs);
764 template<
class AccessorT>
766 correctCellSigns(uint8_t& signs,
768 const AccessorT& acc,
770 const typename AccessorT::ValueType iso)
775 if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 3) signs = uint8_t(~signs);
779 if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 4) signs = uint8_t(~signs);
783 if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 1) signs = uint8_t(~signs);
787 if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 2) signs = uint8_t(~signs);
791 if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 6) signs = uint8_t(~signs);
795 if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 5) signs = uint8_t(~signs);
803 template<
class AccessorT>
805 isNonManifold(
const AccessorT& accessor,
const Coord& ijk,
806 typename AccessorT::ValueType isovalue,
const int dim)
808 const int hDim = dim >> 1;
812 p[0] = isInsideValue(accessor.getValue(coord), isovalue);
814 p[1] = isInsideValue(accessor.getValue(coord), isovalue);
816 p[2] = isInsideValue(accessor.getValue(coord), isovalue);
818 p[3] = isInsideValue(accessor.getValue(coord), isovalue);
819 coord[1] += dim; coord[2] = ijk[2];
820 p[4] = isInsideValue(accessor.getValue(coord), isovalue);
822 p[5] = isInsideValue(accessor.getValue(coord), isovalue);
824 p[6] = isInsideValue(accessor.getValue(coord), isovalue);
826 p[7] = isInsideValue(accessor.getValue(coord), isovalue);
830 if (p[0]) signs |= 1u;
831 if (p[1]) signs |= 2u;
832 if (p[2]) signs |= 4u;
833 if (p[3]) signs |= 8u;
834 if (p[4]) signs |= 16u;
835 if (p[5]) signs |= 32u;
836 if (p[6]) signs |= 64u;
837 if (p[7]) signs |= 128u;
838 if (!sAdaptable[signs])
return true;
843 const int i = ijk[0], ip = ijk[0] + hDim, ipp = ijk[0] + dim;
844 const int j = ijk[1], jp = ijk[1] + hDim, jpp = ijk[1] + dim;
845 const int k = ijk[2], kp = ijk[2] + hDim, kpp = ijk[2] + dim;
848 coord.reset(ip, j, k);
849 m = isInsideValue(accessor.getValue(coord), isovalue);
850 if (p[0] != m && p[1] != m)
return true;
853 coord.reset(ipp, j, kp);
854 m = isInsideValue(accessor.getValue(coord), isovalue);
855 if (p[1] != m && p[2] != m)
return true;
858 coord.reset(ip, j, kpp);
859 m = isInsideValue(accessor.getValue(coord), isovalue);
860 if (p[2] != m && p[3] != m)
return true;
863 coord.reset(i, j, kp);
864 m = isInsideValue(accessor.getValue(coord), isovalue);
865 if (p[0] != m && p[3] != m)
return true;
868 coord.reset(ip, jpp, k);
869 m = isInsideValue(accessor.getValue(coord), isovalue);
870 if (p[4] != m && p[5] != m)
return true;
873 coord.reset(ipp, jpp, kp);
874 m = isInsideValue(accessor.getValue(coord), isovalue);
875 if (p[5] != m && p[6] != m)
return true;
878 coord.reset(ip, jpp, kpp);
879 m = isInsideValue(accessor.getValue(coord), isovalue);
880 if (p[6] != m && p[7] != m)
return true;
883 coord.reset(i, jpp, kp);
884 m = isInsideValue(accessor.getValue(coord), isovalue);
885 if (p[7] != m && p[4] != m)
return true;
888 coord.reset(i, jp, k);
889 m = isInsideValue(accessor.getValue(coord), isovalue);
890 if (p[0] != m && p[4] != m)
return true;
893 coord.reset(ipp, jp, k);
894 m = isInsideValue(accessor.getValue(coord), isovalue);
895 if (p[1] != m && p[5] != m)
return true;
898 coord.reset(ipp, jp, kpp);
899 m = isInsideValue(accessor.getValue(coord), isovalue);
900 if (p[2] != m && p[6] != m)
return true;
904 coord.reset(i, jp, kpp);
905 m = isInsideValue(accessor.getValue(coord), isovalue);
906 if (p[3] != m && p[7] != m)
return true;
912 coord.reset(ip, jp, k);
913 m = isInsideValue(accessor.getValue(coord), isovalue);
914 if (p[0] != m && p[1] != m && p[4] != m && p[5] != m)
return true;
917 coord.reset(ipp, jp, kp);
918 m = isInsideValue(accessor.getValue(coord), isovalue);
919 if (p[1] != m && p[2] != m && p[5] != m && p[6] != m)
return true;
922 coord.reset(ip, jp, kpp);
923 m = isInsideValue(accessor.getValue(coord), isovalue);
924 if (p[2] != m && p[3] != m && p[6] != m && p[7] != m)
return true;
927 coord.reset(i, jp, kp);
928 m = isInsideValue(accessor.getValue(coord), isovalue);
929 if (p[0] != m && p[3] != m && p[4] != m && p[7] != m)
return true;
932 coord.reset(ip, j, kp);
933 m = isInsideValue(accessor.getValue(coord), isovalue);
934 if (p[0] != m && p[1] != m && p[2] != m && p[3] != m)
return true;
937 coord.reset(ip, jpp, kp);
938 m = isInsideValue(accessor.getValue(coord), isovalue);
939 if (p[4] != m && p[5] != m && p[6] != m && p[7] != m)
return true;
942 coord.reset(ip, jp, kp);
943 m = isInsideValue(accessor.getValue(coord), isovalue);
944 if (p[0] != m && p[1] != m && p[2] != m && p[3] != m &&
945 p[4] != m && p[5] != m && p[6] != m && p[7] != m)
return true;
954 template <
class LeafType>
956 mergeVoxels(LeafType& leaf,
const Coord& start,
const int dim,
const int regionId)
959 const Coord end = start.offsetBy(dim);
961 for (ijk[0] = start[0]; ijk[0] < end[0]; ++ijk[0]) {
962 for (ijk[1] = start[1]; ijk[1] < end[1]; ++ijk[1]) {
963 for (ijk[2] = start[2]; ijk[2] < end[2]; ++ijk[2]) {
964 leaf.setValueOnly(ijk, regionId);
973 template <
class LeafType>
975 isMergeable(
const LeafType& leaf,
978 typename LeafType::ValueType::value_type adaptivity)
980 if (adaptivity < 1e-6)
return false;
982 using VecT =
typename LeafType::ValueType;
984 const Coord end = start.offsetBy(dim);
986 std::vector<VecT> norms;
987 for (ijk[0] = start[0]; ijk[0] < end[0]; ++ijk[0]) {
988 for (ijk[1] = start[1]; ijk[1] < end[1]; ++ijk[1]) {
989 for (ijk[2] = start[2]; ijk[2] < end[2]; ++ijk[2]) {
990 if (!leaf.isValueOn(ijk))
continue;
991 norms.push_back(leaf.getValue(ijk));
996 const size_t N = norms.size();
997 for (
size_t ni = 0; ni < N; ++ni) {
998 VecT n_i = norms[ni];
999 for (
size_t nj = 0; nj < N; ++nj) {
1000 VecT n_j = norms[nj];
1001 if ((1.0 - n_i.dot(n_j)) > adaptivity)
return false;
1012 inline double evalZeroCrossing(
double v0,
double v1,
double iso) {
return (iso - v0) / (v1 - v0); }
1017 computePoint(
const std::array<double, 8>& values,
1018 const unsigned char signs,
1019 const unsigned char edgeGroup,
1022 Vec3d avg(0.0, 0.0, 0.0);
1025 if (sEdgeGroupTable[signs][1] == edgeGroup) {
1026 avg[0] += evalZeroCrossing(values[0], values[1], iso);
1030 if (sEdgeGroupTable[signs][2] == edgeGroup) {
1032 avg[2] += evalZeroCrossing(values[1], values[2], iso);
1036 if (sEdgeGroupTable[signs][3] == edgeGroup) {
1037 avg[0] += evalZeroCrossing(values[3], values[2], iso);
1042 if (sEdgeGroupTable[signs][4] == edgeGroup) {
1043 avg[2] += evalZeroCrossing(values[0], values[3], iso);
1047 if (sEdgeGroupTable[signs][5] == edgeGroup) {
1048 avg[0] += evalZeroCrossing(values[4], values[5], iso);
1053 if (sEdgeGroupTable[signs][6] == edgeGroup) {
1056 avg[2] += evalZeroCrossing(values[5], values[6], iso);
1060 if (sEdgeGroupTable[signs][7] == edgeGroup) {
1061 avg[0] += evalZeroCrossing(values[7], values[6], iso);
1067 if (sEdgeGroupTable[signs][8] == edgeGroup) {
1069 avg[2] += evalZeroCrossing(values[4], values[7], iso);
1073 if (sEdgeGroupTable[signs][9] == edgeGroup) {
1074 avg[1] += evalZeroCrossing(values[0], values[4], iso);
1078 if (sEdgeGroupTable[signs][10] == edgeGroup) {
1080 avg[1] += evalZeroCrossing(values[1], values[5], iso);
1084 if (sEdgeGroupTable[signs][11] == edgeGroup) {
1086 avg[1] += evalZeroCrossing(values[2], values[6], iso);
1091 if (sEdgeGroupTable[signs][12] == edgeGroup) {
1092 avg[1] += evalZeroCrossing(values[3], values[7], iso);
1098 const double w = 1.0 / double(samples);
1109 computeMaskedPoint(Vec3d& avg,
1110 const std::array<double, 8>& values,
1111 const unsigned char signs,
1112 const unsigned char signsMask,
1113 const unsigned char edgeGroup,
1116 avg = Vec3d(0.0, 0.0, 0.0);
1119 if (sEdgeGroupTable[signs][1] == edgeGroup &&
1120 sEdgeGroupTable[signsMask][1] == 0) {
1121 avg[0] += evalZeroCrossing(values[0], values[1], iso);
1125 if (sEdgeGroupTable[signs][2] == edgeGroup &&
1126 sEdgeGroupTable[signsMask][2] == 0) {
1128 avg[2] += evalZeroCrossing(values[1], values[2], iso);
1132 if (sEdgeGroupTable[signs][3] == edgeGroup &&
1133 sEdgeGroupTable[signsMask][3] == 0) {
1134 avg[0] += evalZeroCrossing(values[3], values[2], iso);
1139 if (sEdgeGroupTable[signs][4] == edgeGroup &&
1140 sEdgeGroupTable[signsMask][4] == 0) {
1141 avg[2] += evalZeroCrossing(values[0], values[3], iso);
1145 if (sEdgeGroupTable[signs][5] == edgeGroup &&
1146 sEdgeGroupTable[signsMask][5] == 0) {
1147 avg[0] += evalZeroCrossing(values[4], values[5], iso);
1152 if (sEdgeGroupTable[signs][6] == edgeGroup &&
1153 sEdgeGroupTable[signsMask][6] == 0) {
1156 avg[2] += evalZeroCrossing(values[5], values[6], iso);
1160 if (sEdgeGroupTable[signs][7] == edgeGroup &&
1161 sEdgeGroupTable[signsMask][7] == 0) {
1162 avg[0] += evalZeroCrossing(values[7], values[6], iso);
1168 if (sEdgeGroupTable[signs][8] == edgeGroup &&
1169 sEdgeGroupTable[signsMask][8] == 0) {
1171 avg[2] += evalZeroCrossing(values[4], values[7], iso);
1175 if (sEdgeGroupTable[signs][9] == edgeGroup &&
1176 sEdgeGroupTable[signsMask][9] == 0) {
1177 avg[1] += evalZeroCrossing(values[0], values[4], iso);
1181 if (sEdgeGroupTable[signs][10] == edgeGroup &&
1182 sEdgeGroupTable[signsMask][10] == 0) {
1184 avg[1] += evalZeroCrossing(values[1], values[5], iso);
1188 if (sEdgeGroupTable[signs][11] == edgeGroup &&
1189 sEdgeGroupTable[signsMask][11] == 0) {
1191 avg[1] += evalZeroCrossing(values[2], values[6], iso);
1196 if (sEdgeGroupTable[signs][12] == edgeGroup &&
1197 sEdgeGroupTable[signsMask][12] == 0) {
1198 avg[1] += evalZeroCrossing(values[3], values[7], iso);
1204 const double w = 1.0 / double(samples);
1215 computeWeightedPoint(
const Vec3d& p,
1216 const std::array<double, 8>& values,
1217 const unsigned char signs,
1218 const unsigned char edgeGroup,
1221 std::vector<Vec3d> samples;
1224 Vec3d avg(0.0, 0.0, 0.0);
1226 if (sEdgeGroupTable[signs][1] == edgeGroup) {
1227 avg[0] = evalZeroCrossing(values[0], values[1], iso);
1231 samples.push_back(avg);
1234 if (sEdgeGroupTable[signs][2] == edgeGroup) {
1237 avg[2] = evalZeroCrossing(values[1], values[2], iso);
1239 samples.push_back(avg);
1242 if (sEdgeGroupTable[signs][3] == edgeGroup) {
1243 avg[0] = evalZeroCrossing(values[3], values[2], iso);
1247 samples.push_back(avg);
1250 if (sEdgeGroupTable[signs][4] == edgeGroup) {
1253 avg[2] = evalZeroCrossing(values[0], values[3], iso);
1255 samples.push_back(avg);
1258 if (sEdgeGroupTable[signs][5] == edgeGroup) {
1259 avg[0] = evalZeroCrossing(values[4], values[5], iso);
1263 samples.push_back(avg);
1266 if (sEdgeGroupTable[signs][6] == edgeGroup) {
1269 avg[2] = evalZeroCrossing(values[5], values[6], iso);
1271 samples.push_back(avg);
1274 if (sEdgeGroupTable[signs][7] == edgeGroup) {
1275 avg[0] = evalZeroCrossing(values[7], values[6], iso);
1279 samples.push_back(avg);
1282 if (sEdgeGroupTable[signs][8] == edgeGroup) {
1285 avg[2] = evalZeroCrossing(values[4], values[7], iso);
1287 samples.push_back(avg);
1290 if (sEdgeGroupTable[signs][9] == edgeGroup) {
1292 avg[1] = evalZeroCrossing(values[0], values[4], iso);
1295 samples.push_back(avg);
1298 if (sEdgeGroupTable[signs][10] == edgeGroup) {
1300 avg[1] = evalZeroCrossing(values[1], values[5], iso);
1303 samples.push_back(avg);
1306 if (sEdgeGroupTable[signs][11] == edgeGroup) {
1308 avg[1] = evalZeroCrossing(values[2], values[6], iso);
1311 samples.push_back(avg);
1314 if (sEdgeGroupTable[signs][12] == edgeGroup) {
1316 avg[1] = evalZeroCrossing(values[3], values[7], iso);
1319 samples.push_back(avg);
1323 if (samples.size() == 1) {
1324 return samples.front();
1327 std::vector<double> weights;
1328 weights.reserve(samples.size());
1330 for (
const Vec3d& s : samples) {
1331 weights.emplace_back((s-p).lengthSqr());
1334 double minWeight = weights.front();
1335 double maxWeight = weights.front();
1337 for (
size_t i = 1, I = weights.size(); i < I; ++i) {
1338 minWeight =
std::min(minWeight, weights[i]);
1339 maxWeight =
std::max(maxWeight, weights[i]);
1342 const double offset = maxWeight + minWeight * 0.1;
1343 for (
size_t i = 0, I = weights.size(); i < I; ++i) {
1344 weights[i] = offset - weights[i];
1347 double weightSum = 0.0;
1348 for (
size_t i = 0, I = weights.size(); i < I; ++i) {
1349 weightSum += weights[i];
1353 for (
size_t i = 0, I = samples.size(); i < I; ++i) {
1354 avg += samples[i] * (weights[i] / weightSum);
1364 computeCellPoints(std::array<Vec3d, 4>& points,
1365 const std::array<double, 8>& values,
1366 const unsigned char signs,
1370 for (
size_t n = 1, N = sEdgeGroupTable[signs][0] + 1; n < N; ++n, ++offset) {
1372 points[offset] = computePoint(values, signs, uint8_t(n), iso);
1382 matchEdgeGroup(
unsigned char groupId,
unsigned char lhsSigns,
unsigned char rhsSigns)
1385 for (
size_t i = 1; i <= 12; ++i) {
1386 if (sEdgeGroupTable[lhsSigns][i] == groupId && sEdgeGroupTable[rhsSigns][i] != 0) {
1387 id = sEdgeGroupTable[rhsSigns][i];
1400 computeCellPoints(std::array<Vec3d, 4>& points,
1401 std::array<bool, 4>& weightedPointMask,
1402 const std::array<double, 8>& lhsValues,
1403 const std::array<double, 8>& rhsValues,
1404 const unsigned char lhsSigns,
1405 const unsigned char rhsSigns,
1407 const size_t pointIdx,
1408 const uint32_t * seamPointArray)
1411 for (
size_t n = 1, N = sEdgeGroupTable[lhsSigns][0] + 1; n < N; ++n, ++offset)
1414 const int id = matchEdgeGroup(uint8_t(n), lhsSigns, rhsSigns);
1418 const unsigned char e = uint8_t(
id);
1419 const uint32_t quantizedPoint = seamPointArray[pointIdx + (
id - 1)];
1421 if ((quantizedPoint & MASK_DIRTY_BIT) && !(quantizedPoint & MASK_INVALID_BIT)) {
1422 const Vec3d p = unpackPoint(quantizedPoint);
1423 points[offset] = computeWeightedPoint(p, rhsValues, rhsSigns, e, iso);
1424 weightedPointMask[offset] =
true;
1426 points[offset] = computePoint(rhsValues, rhsSigns, e, iso);
1427 weightedPointMask[offset] =
false;
1431 points[offset] = computePoint(lhsValues, lhsSigns, uint8_t(n), iso);
1432 weightedPointMask[offset] =
false;
1439 template <
typename InputTreeType>
1440 struct ComputePoints
1442 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
1443 using InputValueType =
typename InputLeafNodeType::ValueType;
1445 using Int16TreeType =
typename InputTreeType::template ValueConverter<Int16>::Type;
1446 using Int16LeafNodeType =
typename Int16TreeType::LeafNodeType;
1448 using Index32TreeType =
typename InputTreeType::template ValueConverter<Index32>::Type;
1449 using Index32LeafNodeType =
typename Index32TreeType::LeafNodeType;
1451 ComputePoints(
Vec3s * pointArray,
1452 const InputTreeType& inputTree,
1453 const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
1454 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
1455 const std::unique_ptr<
Index32[]>& leafNodeOffsets,
1459 void setRefData(
const InputTreeType& refInputTree,
1460 const Index32TreeType& refPointIndexTree,
1461 const Int16TreeType& refSignFlagsTree,
1462 const uint32_t * quantizedSeamLinePoints,
1463 uint8_t * seamLinePointsFlags);
1465 void operator()(
const tbb::blocked_range<size_t>&)
const;
1468 Vec3s *
const mPoints;
1469 InputTreeType
const *
const mInputTree;
1470 Index32LeafNodeType *
const *
const mPointIndexNodes;
1471 Int16LeafNodeType
const *
const *
const mSignFlagsNodes;
1472 Index32 const *
const mNodeOffsets;
1474 double const mIsovalue;
1476 InputTreeType
const * mRefInputTree;
1477 Index32TreeType
const * mRefPointIndexTree;
1478 Int16TreeType
const * mRefSignFlagsTree;
1479 uint32_t
const * mQuantizedSeamLinePoints;
1480 uint8_t * mSeamLinePointsFlags;
1484 template <
typename InputTreeType>
1485 ComputePoints<InputTreeType>::ComputePoints(
1487 const InputTreeType& inputTree,
1488 const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
1489 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
1490 const std::unique_ptr<
Index32[]>& leafNodeOffsets,
1493 : mPoints(pointArray)
1494 , mInputTree(&inputTree)
1495 , mPointIndexNodes(pointIndexLeafNodes.data())
1496 , mSignFlagsNodes(signFlagsLeafNodes.data())
1497 , mNodeOffsets(leafNodeOffsets.get())
1500 , mRefInputTree(
nullptr)
1501 , mRefPointIndexTree(
nullptr)
1502 , mRefSignFlagsTree(
nullptr)
1503 , mQuantizedSeamLinePoints(
nullptr)
1504 , mSeamLinePointsFlags(
nullptr)
1508 template <
typename InputTreeType>
1510 ComputePoints<InputTreeType>::setRefData(
1511 const InputTreeType& refInputTree,
1512 const Index32TreeType& refPointIndexTree,
1513 const Int16TreeType& refSignFlagsTree,
1514 const uint32_t * quantizedSeamLinePoints,
1515 uint8_t * seamLinePointsFlags)
1517 mRefInputTree = &refInputTree;
1518 mRefPointIndexTree = &refPointIndexTree;
1519 mRefSignFlagsTree = &refSignFlagsTree;
1520 mQuantizedSeamLinePoints = quantizedSeamLinePoints;
1521 mSeamLinePointsFlags = seamLinePointsFlags;
1524 template <
typename InputTreeType>
1526 ComputePoints<InputTreeType>::operator()(
const tbb::blocked_range<size_t>& range)
const 1532 using IndexType =
typename Index32TreeType::ValueType;
1534 using IndexArray = std::vector<Index>;
1535 using IndexArrayMap = std::map<IndexType, IndexArray>;
1537 InputTreeAccessor inputAcc(*mInputTree);
1541 std::array<Vec3d, 4> points;
1542 std::array<bool, 4> weightedPointMask;
1543 std::array<double, 8> values, refValues;
1544 const double iso = mIsovalue;
1548 std::unique_ptr<InputTreeAccessor> refInputAcc;
1549 std::unique_ptr<Index32TreeAccessor> refPointIndexAcc;
1550 std::unique_ptr<Int16TreeAccessor> refSignFlagsAcc;
1552 const bool hasReferenceData = mRefInputTree && mRefPointIndexTree && mRefSignFlagsTree;
1554 if (hasReferenceData) {
1555 refInputAcc.reset(
new InputTreeAccessor(*mRefInputTree));
1556 refPointIndexAcc.reset(
new Index32TreeAccessor(*mRefPointIndexTree));
1557 refSignFlagsAcc.reset(
new Int16TreeAccessor(*mRefSignFlagsTree));
1560 for (
size_t n = range.begin(), N = range.end(); n != N; ++n)
1562 Index32LeafNodeType& pointIndexNode = *mPointIndexNodes[n];
1563 const Coord& origin = pointIndexNode.origin();
1565 const Int16LeafNodeType& signFlagsNode = *mSignFlagsNodes[n];
1566 const InputLeafNodeType * inputNode = inputAcc.probeConstLeaf(origin);
1569 const InputLeafNodeType * refInputNode =
nullptr;
1570 const Index32LeafNodeType * refPointIndexNode =
nullptr;
1571 const Int16LeafNodeType * refSignFlagsNode =
nullptr;
1573 if (hasReferenceData) {
1574 refInputNode = refInputAcc->probeConstLeaf(origin);
1575 refPointIndexNode = refPointIndexAcc->probeConstLeaf(origin);
1576 refSignFlagsNode = refSignFlagsAcc->probeConstLeaf(origin);
1579 IndexType pointOffset = IndexType(mNodeOffsets[n]);
1580 IndexArrayMap regions;
1582 auto*
const pidxData = pointIndexNode.buffer().data();
1583 const auto*
const sfData = signFlagsNode.buffer().data();
1585 for (
auto it = pointIndexNode.beginValueOn(); it; ++it)
1587 const Index offset = it.pos();
1588 IndexType&
id = pidxData[offset];
1592 regions[id].push_back(offset);
1599 const Int16 flags = sfData[offset];
1600 const uint8_t signs = uint8_t(SIGNS & flags);
1601 uint8_t refSigns = 0;
1603 if ((flags & SEAM) && refPointIndexNode && refSignFlagsNode) {
1604 if (refSignFlagsNode->isValueOn(offset)) {
1605 refSigns = uint8_t(SIGNS & refSignFlagsNode->getValue(offset));
1609 ijk = Index32LeafNodeType::offsetToLocalCoord(offset);
1611 const bool inclusiveCell = inputNode && isInternalLeafCoord<InputLeafNodeType>(ijk);
1615 if (inclusiveCell) getCellVertexValues(*inputNode, offset, values);
1616 else getCellVertexValues(inputAcc, ijk, values);
1618 size_t count, weightcount;
1620 if (refSigns == 0) {
1621 count = computeCellPoints(points, values, signs, iso);
1624 if (inclusiveCell && refInputNode) {
1625 getCellVertexValues(*refInputNode, offset, refValues);
1627 getCellVertexValues(*refInputAcc, ijk, refValues);
1629 count = computeCellPoints(points, weightedPointMask, values, refValues, signs, refSigns,
1630 iso, refPointIndexNode->getValue(offset), mQuantizedSeamLinePoints);
1631 weightcount = count;
1634 xyz = ijk.asVec3d();
1636 for (
size_t i = 0; i < count; ++i) {
1638 Vec3d& point = points[i];
1641 if (!std::isfinite(point[0]) ||
1642 !std::isfinite(point[1]) ||
1643 !std::isfinite(point[2]))
1646 "VolumeToMesh encountered NaNs or infs in the input VDB!" 1647 " Hint: Check the input and consider using the \"Diagnostics\" tool " 1648 "to detect and resolve the NaNs.");
1652 point = mTransform.indexToWorld(point);
1654 Vec3s& pos = mPoints[pointOffset];
1655 pos[0] = float(point[0]);
1656 pos[1] = float(point[1]);
1657 pos[2] = float(point[2]);
1659 if (mSeamLinePointsFlags && weightcount && weightedPointMask[i]) {
1660 mSeamLinePointsFlags[pointOffset] = uint8_t(1);
1668 for (
auto it = regions.begin(); it != regions.end(); ++it)
1672 const IndexArray& voxels = it->second;
1673 for (
size_t i = 0, I = voxels.size(); i < I; ++i) {
1675 const Index offset = voxels[i];
1676 ijk = Index32LeafNodeType::offsetToLocalCoord(offset);
1678 const bool inclusiveCell = inputNode && isInternalLeafCoord<InputLeafNodeType>(ijk);
1682 pidxData[offset] = pointOffset;
1684 const uint8_t signs = uint8_t(SIGNS & sfData[offset]);
1686 if (inclusiveCell) getCellVertexValues(*inputNode, offset, values);
1687 else getCellVertexValues(inputAcc, ijk, values);
1689 computeCellPoints(points, values, signs, iso);
1691 avg[0] += double(ijk[0]) + points[0][0];
1692 avg[1] += double(ijk[1]) + points[0][1];
1693 avg[2] += double(ijk[2]) + points[0][2];
1696 if (voxels.size() > 1) {
1697 const double w = 1.0 / double(voxels.size());
1701 avg = mTransform.indexToWorld(avg);
1703 Vec3s& pos = mPoints[pointOffset];
1704 pos[0] = float(avg[0]);
1705 pos[1] = float(avg[1]);
1706 pos[2] = float(avg[2]);
1717 template <
typename InputTreeType>
1718 struct SeamLineWeights
1720 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
1721 using InputValueType =
typename InputLeafNodeType::ValueType;
1723 using Int16TreeType =
typename InputTreeType::template ValueConverter<Int16>::Type;
1724 using Int16LeafNodeType =
typename Int16TreeType::LeafNodeType;
1726 using Index32TreeType =
typename InputTreeType::template ValueConverter<Index32>::Type;
1727 using Index32LeafNodeType =
typename Index32TreeType::LeafNodeType;
1729 SeamLineWeights(
const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
1730 const InputTreeType& inputTree,
1731 const Index32TreeType& refPointIndexTree,
1732 const Int16TreeType& refSignFlagsTree,
1733 uint32_t * quantizedPoints,
1735 : mSignFlagsNodes(signFlagsLeafNodes.data())
1736 , mInputTree(&inputTree)
1737 , mRefPointIndexTree(&refPointIndexTree)
1738 , mRefSignFlagsTree(&refSignFlagsTree)
1739 , mQuantizedPoints(quantizedPoints)
1744 void operator()(
const tbb::blocked_range<size_t>& range)
const 1750 std::array<double, 8> values;
1751 const double iso = double(mIsovalue);
1755 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
1757 const Int16LeafNodeType& signFlagsNode = *mSignFlagsNodes[n];
1758 const Coord& origin = signFlagsNode.origin();
1760 const Int16LeafNodeType * refSignNode = signFlagsTreeAcc.
probeConstLeaf(origin);
1761 if (!refSignNode)
continue;
1763 const Index32LeafNodeType* refPointIndexNode =
1765 if (!refPointIndexNode)
continue;
1767 const InputLeafNodeType* inputNode = inputTreeAcc.
probeConstLeaf(origin);
1769 const auto*
const sfData = signFlagsNode.buffer().data();
1770 const auto*
const rfIdxData = refPointIndexNode->buffer().data();
1771 const auto*
const rsData = refSignNode->buffer().data();
1773 for (
auto it = signFlagsNode.cbeginValueOn(); it; ++it)
1775 const Index offset = it.pos();
1776 const Int16 flags = sfData[offset];
1778 ijk = Index32LeafNodeType::offsetToLocalCoord(offset);
1780 const bool inclusiveCell = inputNode && isInternalLeafCoord<InputLeafNodeType>(ijk);
1784 if ((flags & SEAM) && refSignNode->isValueOn(offset)) {
1786 const uint8_t lhsSigns = uint8_t(SIGNS & flags);
1787 const uint8_t rhsSigns = uint8_t(SIGNS & rsData[offset]);
1789 if (inclusiveCell) getCellVertexValues(*inputNode, offset, values);
1790 else getCellVertexValues(inputTreeAcc, ijk, values);
1792 for (
unsigned i = 1, I = sEdgeGroupTable[lhsSigns][0] + 1; i < I; ++i) {
1794 const int id = matchEdgeGroup(uint8_t(i), lhsSigns, rhsSigns);
1798 uint32_t& data = mQuantizedPoints[rfIdxData[offset] + (
id - 1)];
1800 if (!(data & MASK_DIRTY_BIT)) {
1802 const int samples = computeMaskedPoint(
1803 pos, values, lhsSigns, rhsSigns, uint8_t(i), iso);
1805 if (samples > 0) data = packPoint(pos);
1806 else data = MASK_INVALID_BIT;
1808 data |= MASK_DIRTY_BIT;
1818 Int16LeafNodeType
const *
const *
const mSignFlagsNodes;
1819 InputTreeType
const *
const mInputTree;
1820 Index32TreeType
const *
const mRefPointIndexTree;
1821 Int16TreeType
const *
const mRefSignFlagsTree;
1822 uint32_t *
const mQuantizedPoints;
1823 InputValueType
const mIsovalue;
1827 template <
typename TreeType>
1828 struct SetSeamLineFlags
1830 using LeafNodeType =
typename TreeType::LeafNodeType;
1832 SetSeamLineFlags(
const std::vector<LeafNodeType*>& signFlagsLeafNodes,
1833 const TreeType& refSignFlagsTree)
1834 : mSignFlagsNodes(signFlagsLeafNodes.data())
1835 , mRefSignFlagsTree(&refSignFlagsTree)
1839 void operator()(
const tbb::blocked_range<size_t>& range)
const 1843 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
1845 LeafNodeType& signFlagsNode = *mSignFlagsNodes[n];
1846 const Coord& origin = signFlagsNode.origin();
1848 const LeafNodeType* refSignNode = refSignFlagsTreeAcc.
probeConstLeaf(origin);
1849 if (!refSignNode)
continue;
1851 const auto*
const rsData = refSignNode->buffer().data();
1852 auto*
const sfData = signFlagsNode.buffer().data();
1854 for (
auto it = signFlagsNode.cbeginValueOn(); it; ++it) {
1855 const Index offset = it.pos();
1857 const uint8_t rhsSigns = uint8_t(rsData[offset] & SIGNS);
1859 if (sEdgeGroupTable[rhsSigns][0] > 0) {
1861 typename LeafNodeType::ValueType& value = sfData[offset];
1862 const uint8_t lhsSigns = uint8_t(value & SIGNS);
1864 if (rhsSigns != lhsSigns) {
1875 LeafNodeType *
const *
const mSignFlagsNodes;
1876 TreeType
const *
const mRefSignFlagsTree;
1880 template <
typename BoolTreeType,
typename SignDataType>
1881 struct TransferSeamLineFlags
1883 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
1885 using SignDataTreeType =
typename BoolTreeType::template ValueConverter<SignDataType>::Type;
1886 using SignDataLeafNodeType =
typename SignDataTreeType::LeafNodeType;
1888 TransferSeamLineFlags(
const std::vector<SignDataLeafNodeType*>& signFlagsLeafNodes,
1889 const BoolTreeType& maskTree)
1890 : mSignFlagsNodes(signFlagsLeafNodes.data())
1891 , mMaskTree(&maskTree)
1895 void operator()(
const tbb::blocked_range<size_t>& range)
const 1899 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
1901 SignDataLeafNodeType& signFlagsNode = *mSignFlagsNodes[n];
1902 const Coord& origin = signFlagsNode.origin();
1904 const BoolLeafNodeType * maskNode = maskAcc.
probeConstLeaf(origin);
1905 if (!maskNode)
continue;
1907 auto*
const sfData = signFlagsNode.buffer().data();
1909 for (
auto it = signFlagsNode.cbeginValueOn(); it; ++it) {
1910 const Index offset = it.pos();
1912 if (maskNode->isValueOn(offset)) {
1913 sfData[offset] |= SEAM;
1920 SignDataLeafNodeType *
const *
const mSignFlagsNodes;
1921 BoolTreeType
const *
const mMaskTree;
1925 template <
typename TreeType>
1926 struct MaskSeamLineVoxels
1928 using LeafNodeType =
typename TreeType::LeafNodeType;
1929 using BoolTreeType =
typename TreeType::template ValueConverter<bool>::Type;
1931 MaskSeamLineVoxels(
const std::vector<LeafNodeType*>& signFlagsLeafNodes,
1932 const TreeType& signFlagsTree,
1934 : mSignFlagsNodes(signFlagsLeafNodes.data())
1935 , mSignFlagsTree(&signFlagsTree)
1941 MaskSeamLineVoxels(MaskSeamLineVoxels& rhs,
tbb::split)
1942 : mSignFlagsNodes(rhs.mSignFlagsNodes)
1943 , mSignFlagsTree(rhs.mSignFlagsTree)
1949 void join(MaskSeamLineVoxels& rhs) { mMask->merge(*rhs.mMask); }
1951 void operator()(
const tbb::blocked_range<size_t>& range)
1953 using ValueOnCIter =
typename LeafNodeType::ValueOnCIter;
1954 using ValueType =
typename LeafNodeType::ValueType;
1960 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
1962 LeafNodeType& signFlagsNode = *mSignFlagsNodes[n];
1963 auto*
const sfData = signFlagsNode.buffer().data();
1965 for (ValueOnCIter it = signFlagsNode.cbeginValueOn(); it; ++it) {
1967 const ValueType flags = sfData[it.pos()];
1969 if (!(flags & SEAM) && (flags & EDGES)) {
1971 ijk = it.getCoord();
1973 bool isSeamLineVoxel =
false;
1975 if (flags & XEDGE) {
1977 isSeamLineVoxel = (signFlagsAcc.
getValue(ijk) & SEAM);
1979 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) & SEAM);
1981 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) & SEAM);
1985 if (!isSeamLineVoxel && flags & YEDGE) {
1987 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) & SEAM);
1989 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) & SEAM);
1991 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) & SEAM);
1995 if (!isSeamLineVoxel && flags & ZEDGE) {
1997 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) & SEAM);
1999 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) & SEAM);
2001 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) & SEAM);
2005 if (isSeamLineVoxel) {
2006 maskAcc.
setValue(it.getCoord(),
true);
2015 LeafNodeType *
const *
const mSignFlagsNodes;
2016 TreeType
const *
const mSignFlagsTree;
2017 BoolTreeType mTempMask;
2018 BoolTreeType *
const mMask;
2022 template<
typename SignDataTreeType>
2024 markSeamLineData(SignDataTreeType& signFlagsTree,
const SignDataTreeType& refSignFlagsTree)
2026 using SignDataType =
typename SignDataTreeType::ValueType;
2027 using SignDataLeafNodeType =
typename SignDataTreeType::LeafNodeType;
2028 using BoolTreeType =
typename SignDataTreeType::template ValueConverter<bool>::Type;
2030 std::vector<SignDataLeafNodeType*> signFlagsLeafNodes;
2031 signFlagsTree.getNodes(signFlagsLeafNodes);
2033 const tbb::blocked_range<size_t> nodeRange(0, signFlagsLeafNodes.size());
2035 tbb::parallel_for(nodeRange,
2036 SetSeamLineFlags<SignDataTreeType>(signFlagsLeafNodes, refSignFlagsTree));
2038 BoolTreeType seamLineMaskTree(
false);
2040 MaskSeamLineVoxels<SignDataTreeType>
2041 maskSeamLine(signFlagsLeafNodes, signFlagsTree, seamLineMaskTree);
2043 tbb::parallel_reduce(nodeRange, maskSeamLine);
2045 tbb::parallel_for(nodeRange,
2046 TransferSeamLineFlags<BoolTreeType, SignDataType>(signFlagsLeafNodes, seamLineMaskTree));
2053 template <
typename InputGr
idType>
2054 struct MergeVoxelRegions
2056 using InputTreeType =
typename InputGridType::TreeType;
2057 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
2058 using InputValueType =
typename InputLeafNodeType::ValueType;
2060 using FloatTreeType =
typename InputTreeType::template ValueConverter<float>::Type;
2061 using FloatLeafNodeType =
typename FloatTreeType::LeafNodeType;
2064 using Int16TreeType =
typename InputTreeType::template ValueConverter<Int16>::Type;
2065 using Int16LeafNodeType =
typename Int16TreeType::LeafNodeType;
2067 using Index32TreeType =
typename InputTreeType::template ValueConverter<Index32>::Type;
2068 using Index32LeafNodeType =
typename Index32TreeType::LeafNodeType;
2070 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
2071 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
2073 using MaskTreeType =
typename InputTreeType::template ValueConverter<ValueMask>::Type;
2074 using MaskLeafNodeType =
typename MaskTreeType::LeafNodeType;
2076 MergeVoxelRegions(
const InputGridType& inputGrid,
2077 const Index32TreeType& pointIndexTree,
2078 const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
2079 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
2082 bool invertSurfaceOrientation);
2084 void setSpatialAdaptivity(
const FloatGridType& grid)
2086 mSpatialAdaptivityTree = &grid.tree();
2087 mSpatialAdaptivityTransform = &grid.transform();
2090 void setAdaptivityMask(
const BoolTreeType& mask)
2095 void setRefSignFlagsData(
const Int16TreeType& signFlagsData,
float internalAdaptivity)
2097 mRefSignFlagsTree = &signFlagsData;
2098 mInternalAdaptivity = internalAdaptivity;
2101 void operator()(
const tbb::blocked_range<size_t>&)
const;
2104 InputTreeType
const *
const mInputTree;
2107 Index32TreeType
const *
const mPointIndexTree;
2108 Index32LeafNodeType *
const *
const mPointIndexNodes;
2109 Int16LeafNodeType
const *
const *
const mSignFlagsNodes;
2111 InputValueType mIsovalue;
2112 float mSurfaceAdaptivity, mInternalAdaptivity;
2113 bool mInvertSurfaceOrientation;
2115 FloatTreeType
const * mSpatialAdaptivityTree;
2116 BoolTreeType
const * mMaskTree;
2117 Int16TreeType
const * mRefSignFlagsTree;
2122 template <
typename InputGr
idType>
2123 MergeVoxelRegions<InputGridType>::MergeVoxelRegions(
2124 const InputGridType& inputGrid,
2125 const Index32TreeType& pointIndexTree,
2126 const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
2127 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
2130 bool invertSurfaceOrientation)
2131 : mInputTree(&inputGrid.tree())
2132 , mInputTransform(&inputGrid.transform())
2133 , mPointIndexTree(&pointIndexTree)
2134 , mPointIndexNodes(pointIndexLeafNodes.data())
2135 , mSignFlagsNodes(signFlagsLeafNodes.data())
2137 , mSurfaceAdaptivity(adaptivity)
2138 , mInternalAdaptivity(adaptivity)
2139 , mInvertSurfaceOrientation(invertSurfaceOrientation)
2140 , mSpatialAdaptivityTree(
nullptr)
2141 , mMaskTree(
nullptr)
2142 , mRefSignFlagsTree(
nullptr)
2143 , mSpatialAdaptivityTransform(
nullptr)
2148 template <
typename InputGr
idType>
2150 MergeVoxelRegions<InputGridType>::operator()(
const tbb::blocked_range<size_t>& range)
const 2153 using Vec3sLeafNodeType =
typename InputLeafNodeType::template ValueConverter<Vec3sType>::Type;
2161 std::unique_ptr<FloatTreeAccessor> spatialAdaptivityAcc;
2162 if (mSpatialAdaptivityTree && mSpatialAdaptivityTransform) {
2163 spatialAdaptivityAcc.reset(
new FloatTreeAccessor(*mSpatialAdaptivityTree));
2166 std::unique_ptr<BoolTreeAccessor> maskAcc;
2168 maskAcc.reset(
new BoolTreeAccessor(*mMaskTree));
2171 std::unique_ptr<Int16TreeAccessor> refSignFlagsAcc;
2172 if (mRefSignFlagsTree) {
2173 refSignFlagsAcc.reset(
new Int16TreeAccessor(*mRefSignFlagsTree));
2176 InputTreeAccessor inputAcc(*mInputTree);
2177 Index32TreeAccessor pointIndexAcc(*mPointIndexTree);
2179 MaskLeafNodeType mask;
2181 const bool invertGradientDir = mInvertSurfaceOrientation || isBoolValue<InputValueType>();
2182 std::unique_ptr<Vec3sLeafNodeType> gradientNode;
2185 const int LeafDim = InputLeafNodeType::DIM;
2187 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
2189 mask.setValuesOff();
2191 const Int16LeafNodeType& signFlagsNode = *mSignFlagsNodes[n];
2192 Index32LeafNodeType& pointIndexNode = *mPointIndexNodes[n];
2194 const Coord& origin = pointIndexNode.origin();
2195 end = origin.offsetBy(LeafDim);
2199 const BoolLeafNodeType* maskLeaf = maskAcc->probeConstLeaf(origin);
2200 if (maskLeaf !=
nullptr) {
2201 for (
auto it = maskLeaf->cbeginValueOn(); it; ++it)
2203 mask.setActiveState(it.getCoord() & ~1u,
true);
2208 float adaptivity = (refSignFlagsAcc && !refSignFlagsAcc->probeConstLeaf(origin)) ?
2209 mInternalAdaptivity : mSurfaceAdaptivity;
2211 bool useGradients = adaptivity < 1.0f;
2214 FloatLeafNodeType adaptivityLeaf(origin, adaptivity);
2216 if (spatialAdaptivityAcc) {
2217 useGradients =
false;
2218 for (
Index offset = 0; offset < FloatLeafNodeType::NUM_VALUES; ++offset) {
2219 ijk = adaptivityLeaf.offsetToGlobalCoord(offset);
2220 ijk = mSpatialAdaptivityTransform->worldToIndexCellCentered(
2221 mInputTransform->indexToWorld(ijk));
2222 float weight = spatialAdaptivityAcc->getValue(ijk);
2223 float adaptivityValue = weight * adaptivity;
2224 if (adaptivityValue < 1.0f) useGradients =
true;
2225 adaptivityLeaf.setValueOnly(offset, adaptivityValue);
2230 for (
auto it = signFlagsNode.cbeginValueOn(); it; ++it) {
2231 const Int16 flags = it.getValue();
2232 const unsigned char signs =
static_cast<unsigned char>(SIGNS & int(flags));
2234 if ((flags & SEAM) || !sAdaptable[signs] || sEdgeGroupTable[signs][0] > 1) {
2236 mask.setActiveState(it.getCoord() & ~1u,
true);
2238 }
else if (flags & EDGES) {
2240 bool maskRegion =
false;
2242 ijk = it.getCoord();
2243 if (!pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2245 if (!maskRegion && flags & XEDGE) {
2247 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2249 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2251 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2255 if (!maskRegion && flags & YEDGE) {
2257 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2259 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2261 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2265 if (!maskRegion && flags & ZEDGE) {
2267 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2269 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2271 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2276 mask.setActiveState(it.getCoord() & ~1u,
true);
2283 for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
2284 for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
2285 for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
2286 if (!mask.isValueOn(ijk) && isNonManifold(inputAcc, ijk, mIsovalue, dim)) {
2287 mask.setActiveState(ijk,
true);
2298 gradientNode->setValuesOff();
2300 gradientNode.reset(
new Vec3sLeafNodeType());
2303 for (
auto it = signFlagsNode.cbeginValueOn(); it; ++it)
2305 ijk = it.getCoord();
2306 if (!mask.isValueOn(ijk & ~1u))
2311 if (invertGradientDir) {
2315 gradientNode->setValueOn(it.pos(), dir);
2322 for ( ; dim <= LeafDim; dim = dim << 1) {
2323 const unsigned coordMask = ~((dim << 1) - 1);
2324 for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
2325 for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
2326 for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
2328 adaptivity = adaptivityLeaf.getValue(ijk);
2330 if (mask.isValueOn(ijk)
2331 || isNonManifold(inputAcc, ijk, mIsovalue, dim)
2332 || (useGradients && !isMergeable(*gradientNode, ijk, dim, adaptivity)))
2334 mask.setActiveState(ijk & coordMask,
true);
2336 mergeVoxels(pointIndexNode, ijk, dim, regionId++);
2350 struct UniformPrimBuilder
2352 UniformPrimBuilder(): mIdx(0), mPolygonPool(
nullptr) {}
2354 void init(
const size_t upperBound,
PolygonPool& quadPool)
2356 mPolygonPool = &quadPool;
2361 template<
typename IndexType>
2365 mPolygonPool->quad(mIdx) = verts;
2367 Vec4I& quad = mPolygonPool->quad(mIdx);
2373 mPolygonPool->quadFlags(mIdx) = flags;
2379 mPolygonPool->trimQuads(mIdx);
2389 struct AdaptivePrimBuilder
2391 AdaptivePrimBuilder() : mQuadIdx(0), mTriangleIdx(0), mPolygonPool(
nullptr) {}
2393 void init(
const size_t upperBound,
PolygonPool& polygonPool)
2395 mPolygonPool = &polygonPool;
2397 mPolygonPool->resetTriangles(upperBound);
2403 template<
typename IndexType>
2406 if (verts[0] != verts[1] && verts[0] != verts[2] && verts[0] != verts[3]
2407 && verts[1] != verts[2] && verts[1] != verts[3] && verts[2] != verts[3]) {
2408 mPolygonPool->quadFlags(mQuadIdx) = flags;
2409 addQuad(verts, reverse);
2411 verts[0] == verts[3] &&
2412 verts[1] != verts[2] &&
2413 verts[1] != verts[0] &&
2414 verts[2] != verts[0]) {
2415 mPolygonPool->triangleFlags(mTriangleIdx) = flags;
2416 addTriangle(verts[0], verts[1], verts[2], reverse);
2418 verts[1] == verts[2] &&
2419 verts[0] != verts[3] &&
2420 verts[0] != verts[1] &&
2421 verts[3] != verts[1]) {
2422 mPolygonPool->triangleFlags(mTriangleIdx) = flags;
2423 addTriangle(verts[0], verts[1], verts[3], reverse);
2425 verts[0] == verts[1] &&
2426 verts[2] != verts[3] &&
2427 verts[2] != verts[0] &&
2428 verts[3] != verts[0]) {
2429 mPolygonPool->triangleFlags(mTriangleIdx) = flags;
2430 addTriangle(verts[0], verts[2], verts[3], reverse);
2432 verts[2] == verts[3] &&
2433 verts[0] != verts[1] &&
2434 verts[0] != verts[2] &&
2435 verts[1] != verts[2]) {
2436 mPolygonPool->triangleFlags(mTriangleIdx) = flags;
2437 addTriangle(verts[0], verts[1], verts[2], reverse);
2444 mPolygonPool->trimQuads(mQuadIdx,
true);
2445 mPolygonPool->trimTrinagles(mTriangleIdx,
true);
2450 template<
typename IndexType>
2454 mPolygonPool->quad(mQuadIdx) = verts;
2456 Vec4I& quad = mPolygonPool->quad(mQuadIdx);
2465 void addTriangle(
unsigned v0,
unsigned v1,
unsigned v2,
bool reverse)
2467 Vec3I& prim = mPolygonPool->triangle(mTriangleIdx);
2481 size_t mQuadIdx, mTriangleIdx;
2486 template<
typename SignAccT,
typename IdxAccT,
typename PrimBuilder>
2489 bool invertSurfaceOrientation,
2492 const Vec3i& offsets,
2494 const SignAccT& signAcc,
2495 const IdxAccT& idxAcc,
2496 PrimBuilder& mesher)
2498 using IndexType =
typename IdxAccT::ValueType;
2501 const bool isActive = idxAcc.probeValue(ijk, v0);
2508 bool isInside = flags & INSIDE;
2510 isInside = invertSurfaceOrientation ? !isInside : isInside;
2515 if (flags & XEDGE) {
2517 quad[0] = v0 + offsets[0];
2521 bool activeValues = idxAcc.probeValue(coord, quad[1]);
2522 uint8_t cell = uint8_t(SIGNS & signAcc.getValue(coord));
2523 quad[1] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][5] - 1 : 0;
2527 activeValues = activeValues && idxAcc.probeValue(coord, quad[2]);
2528 cell = uint8_t(SIGNS & signAcc.getValue(coord));
2529 quad[2] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][7] - 1 : 0;
2533 activeValues = activeValues && idxAcc.probeValue(coord, quad[3]);
2534 cell = uint8_t(SIGNS & signAcc.getValue(coord));
2535 quad[3] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][3] - 1 : 0;
2538 mesher.addPrim(quad, isInside, tag[
bool(refFlags & XEDGE)]);
2545 if (flags & YEDGE) {
2547 quad[0] = v0 + offsets[1];
2551 bool activeValues = idxAcc.probeValue(coord, quad[1]);
2552 uint8_t cell = uint8_t(SIGNS & signAcc.getValue(coord));
2553 quad[1] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][12] - 1 : 0;
2557 activeValues = activeValues && idxAcc.probeValue(coord, quad[2]);
2558 cell = uint8_t(SIGNS & signAcc.getValue(coord));
2559 quad[2] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][11] - 1 : 0;
2563 activeValues = activeValues && idxAcc.probeValue(coord, quad[3]);
2564 cell = uint8_t(SIGNS & signAcc.getValue(coord));
2565 quad[3] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][10] - 1 : 0;
2568 mesher.addPrim(quad, isInside, tag[
bool(refFlags & YEDGE)]);
2575 if (flags & ZEDGE) {
2577 quad[0] = v0 + offsets[2];
2581 bool activeValues = idxAcc.probeValue(coord, quad[1]);
2582 uint8_t cell = uint8_t(SIGNS & signAcc.getValue(coord));
2583 quad[1] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][8] - 1 : 0;
2587 activeValues = activeValues && idxAcc.probeValue(coord, quad[2]);
2588 cell = uint8_t(SIGNS & signAcc.getValue(coord));
2589 quad[2] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][6] - 1 : 0;
2593 activeValues = activeValues && idxAcc.probeValue(coord, quad[3]);
2594 cell = uint8_t(SIGNS & signAcc.getValue(coord));
2595 quad[3] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][2] - 1 : 0;
2598 mesher.addPrim(quad, !isInside, tag[
bool(refFlags & ZEDGE)]);
2607 template<
typename InputTreeType>
2608 struct MaskTileBorders
2610 using InputValueType =
typename InputTreeType::ValueType;
2611 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
2614 MaskTileBorders(
const InputTreeType& inputTree, InputValueType iso,
2615 BoolTreeType& mask,
const Vec4i* tileArray)
2616 : mInputTree(&inputTree)
2620 , mTileArray(tileArray)
2624 MaskTileBorders(MaskTileBorders& rhs,
tbb::split)
2625 : mInputTree(rhs.mInputTree)
2626 , mIsovalue(rhs.mIsovalue)
2629 , mTileArray(rhs.mTileArray)
2633 void join(MaskTileBorders& rhs) { mMask->merge(*rhs.mMask); }
2635 void operator()(
const tbb::blocked_range<size_t>&);
2638 InputTreeType
const *
const mInputTree;
2639 InputValueType
const mIsovalue;
2640 BoolTreeType mTempMask;
2641 BoolTreeType *
const mMask;
2642 Vec4i const *
const mTileArray;
2646 template<
typename InputTreeType>
2648 MaskTileBorders<InputTreeType>::operator()(
const tbb::blocked_range<size_t>& range)
2652 CoordBBox region, bbox;
2655 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
2657 const Vec4i& tile = mTileArray[n];
2659 bbox.min()[0] = tile[0];
2660 bbox.min()[1] = tile[1];
2661 bbox.min()[2] = tile[2];
2662 bbox.max() = bbox.min();
2663 bbox.max().offset(tile[3]);
2665 InputValueType value = mInputTree->background();
2667 const bool isInside = isInsideValue(inputTreeAcc.
getValue(bbox.min()), mIsovalue);
2668 const int valueDepth = inputTreeAcc.
getValueDepth(bbox.min());
2676 bool processRegion =
true;
2678 processRegion = isInside != isInsideValue(inputTreeAcc.
getValue(nijk), mIsovalue);
2681 if (processRegion) {
2684 region.min()[0] = region.max()[0] = ijk[0];
2685 mMask->fill(region,
false);
2692 processRegion =
true;
2694 processRegion = (!inputTreeAcc.
probeValue(ijk, value)
2695 && isInside != isInsideValue(value, mIsovalue));
2698 if (processRegion) {
2701 region.min()[0] = region.max()[0] = ijk[0];
2702 mMask->fill(region,
false);
2712 processRegion =
true;
2714 processRegion = isInside != isInsideValue(inputTreeAcc.
getValue(nijk), mIsovalue);
2717 if (processRegion) {
2720 region.min()[1] = region.max()[1] = ijk[1];
2721 mMask->fill(region,
false);
2728 processRegion =
true;
2730 processRegion = (!inputTreeAcc.
probeValue(ijk, value)
2731 && isInside != isInsideValue(value, mIsovalue));
2734 if (processRegion) {
2737 region.min()[1] = region.max()[1] = ijk[1];
2738 mMask->fill(region,
false);
2748 processRegion =
true;
2750 processRegion = isInside != isInsideValue(inputTreeAcc.
getValue(nijk), mIsovalue);
2753 if (processRegion) {
2756 region.min()[2] = region.max()[2] = ijk[2];
2757 mMask->fill(region,
false);
2763 processRegion =
true;
2765 processRegion = (!inputTreeAcc.
probeValue(ijk, value)
2766 && isInside != isInsideValue(value, mIsovalue));
2769 if (processRegion) {
2772 region.min()[2] = region.max()[2] = ijk[2];
2773 mMask->fill(region,
false);
2779 template<
typename InputTreeType>
2781 maskActiveTileBorders(
const InputTreeType& inputTree,
2782 const typename InputTreeType::ValueType iso,
2783 typename InputTreeType::template ValueConverter<bool>::Type& mask)
2785 typename InputTreeType::ValueOnCIter tileIter(inputTree);
2786 tileIter.setMaxDepth(InputTreeType::ValueOnCIter::LEAF_DEPTH - 1);
2788 size_t tileCount = 0;
2789 for ( ; tileIter; ++tileIter) {
2793 if (tileCount > 0) {
2794 std::unique_ptr<Vec4i[]> tiles(
new Vec4i[tileCount]);
2799 tileIter = inputTree.cbeginValueOn();
2800 tileIter.setMaxDepth(InputTreeType::ValueOnCIter::LEAF_DEPTH - 1);
2802 for (; tileIter; ++tileIter) {
2803 Vec4i& tile = tiles[index++];
2804 tileIter.getBoundingBox(bbox);
2805 tile[0] = bbox.min()[0];
2806 tile[1] = bbox.min()[1];
2807 tile[2] = bbox.min()[2];
2808 tile[3] = bbox.max()[0] - bbox.min()[0];
2811 MaskTileBorders<InputTreeType>
op(inputTree, iso, mask, tiles.get());
2812 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, tileCount),
op);
2824 PointListCopy(
const PointList& pointsIn, std::vector<Vec3s>& pointsOut)
2825 : mPointsIn(pointsIn) , mPointsOut(pointsOut)
2829 void operator()(
const tbb::blocked_range<size_t>& range)
const 2831 for (
size_t n = range.begin(); n < range.end(); ++n) {
2832 mPointsOut[n] = mPointsIn[n];
2838 std::vector<Vec3s>& mPointsOut;
2846 struct LeafNodeVoxelOffsets
2848 using IndexVector = std::vector<Index>;
2850 template<
typename LeafNodeType>
2851 void constructOffsetList();
2854 const IndexVector& core()
const {
return mCore; }
2858 const IndexVector& minX()
const {
return mMinX; }
2861 const IndexVector& maxX()
const {
return mMaxX; }
2865 const IndexVector& minY()
const {
return mMinY; }
2868 const IndexVector& maxY()
const {
return mMaxY; }
2872 const IndexVector& minZ()
const {
return mMinZ; }
2875 const IndexVector& maxZ()
const {
return mMaxZ; }
2879 const IndexVector& internalNeighborsX()
const {
return mInternalNeighborsX; }
2882 const IndexVector& internalNeighborsY()
const {
return mInternalNeighborsY; }
2885 const IndexVector& internalNeighborsZ()
const {
return mInternalNeighborsZ; }
2889 IndexVector mCore, mMinX, mMaxX, mMinY, mMaxY, mMinZ, mMaxZ,
2890 mInternalNeighborsX, mInternalNeighborsY, mInternalNeighborsZ;
2894 template<
typename LeafNodeType>
2896 LeafNodeVoxelOffsets::constructOffsetList()
2900 mCore.reserve((LeafNodeType::DIM - 2) * (LeafNodeType::DIM - 2));
2902 for (
Index x = 1; x < (LeafNodeType::DIM - 1); ++x) {
2903 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2904 for (
Index y = 1; y < (LeafNodeType::DIM - 1); ++y) {
2905 const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
2906 for (
Index z = 1; z < (LeafNodeType::DIM - 1); ++z) {
2907 mCore.push_back(offsetXY + z);
2913 mInternalNeighborsX.clear();
2914 mInternalNeighborsX.reserve(LeafNodeType::SIZE - (LeafNodeType::DIM * LeafNodeType::DIM));
2916 for (
Index x = 0; x < (LeafNodeType::DIM - 1); ++x) {
2917 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2918 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
2919 const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
2920 for (
Index z = 0; z < LeafNodeType::DIM; ++z) {
2921 mInternalNeighborsX.push_back(offsetXY + z);
2927 mInternalNeighborsY.clear();
2928 mInternalNeighborsY.reserve(LeafNodeType::SIZE - (LeafNodeType::DIM * LeafNodeType::DIM));
2930 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
2931 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2932 for (
Index y = 0; y < (LeafNodeType::DIM - 1); ++y) {
2933 const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
2934 for (
Index z = 0; z < LeafNodeType::DIM; ++z) {
2935 mInternalNeighborsY.push_back(offsetXY + z);
2941 mInternalNeighborsZ.clear();
2942 mInternalNeighborsZ.reserve(LeafNodeType::SIZE - (LeafNodeType::DIM * LeafNodeType::DIM));
2944 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
2945 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2946 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
2947 const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
2948 for (
Index z = 0; z < (LeafNodeType::DIM - 1); ++z) {
2949 mInternalNeighborsZ.push_back(offsetXY + z);
2956 mMinX.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
2958 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
2959 const Index offsetXY = (y << LeafNodeType::LOG2DIM);
2960 for (
Index z = 0; z < LeafNodeType::DIM; ++z) {
2961 mMinX.push_back(offsetXY + z);
2968 mMaxX.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
2970 const Index offsetX = (LeafNodeType::DIM - 1) << (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; ++z) {
2974 mMaxX.push_back(offsetXY + z);
2981 mMinY.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
2983 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
2984 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2985 for (
Index z = 0; z < (LeafNodeType::DIM - 1); ++z) {
2986 mMinY.push_back(offsetX + z);
2993 mMaxY.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
2995 const Index offsetY = (LeafNodeType::DIM - 1) << LeafNodeType::LOG2DIM;
2996 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
2997 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2998 for (
Index z = 0; z < (LeafNodeType::DIM - 1); ++z) {
2999 mMaxY.push_back(offsetX + offsetY + z);
3006 mMinZ.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
3008 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
3009 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
3010 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
3011 const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
3012 mMinZ.push_back(offsetXY);
3019 mMaxZ.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
3021 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
3022 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
3023 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
3024 const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
3025 mMaxZ.push_back(offsetXY + (LeafNodeType::DIM - 1));
3036 template<
typename AccessorT,
int _AXIS>
3037 struct VoxelEdgeAccessor
3039 enum { AXIS = _AXIS };
3042 VoxelEdgeAccessor(AccessorT& _acc) : acc(_acc) {}
3044 void set(Coord ijk) {
3046 acc.setActiveState(ijk);
3048 acc.setActiveState(ijk);
3050 acc.setActiveState(ijk);
3052 acc.setActiveState(ijk);
3053 }
else if (_AXIS == 1) {
3054 acc.setActiveState(ijk);
3056 acc.setActiveState(ijk);
3058 acc.setActiveState(ijk);
3060 acc.setActiveState(ijk);
3062 acc.setActiveState(ijk);
3064 acc.setActiveState(ijk);
3066 acc.setActiveState(ijk);
3068 acc.setActiveState(ijk);
3077 template<
typename VoxelEdgeAcc,
typename LeafNodeT>
3079 evalInternalVoxelEdges(VoxelEdgeAcc& edgeAcc,
3080 const LeafNodeT& leafnode,
3081 const LeafNodeVoxelOffsets& voxels,
3082 const typename LeafNodeT::ValueType iso)
3085 const std::vector<Index>* offsets = &voxels.internalNeighborsZ();
3087 if (VoxelEdgeAcc::AXIS == 0) {
3088 nvo = LeafNodeT::DIM * LeafNodeT::DIM;
3089 offsets = &voxels.internalNeighborsX();
3090 }
else if (VoxelEdgeAcc::AXIS == 1) {
3091 nvo = LeafNodeT::DIM;
3092 offsets = &voxels.internalNeighborsY();
3095 const LeafBufferAccessor<LeafNodeT> lhsAcc(leafnode);
3097 for (
size_t n = 0, N = offsets->size(); n < N; ++n) {
3098 const Index& pos = (*offsets)[n];
3099 const bool isActive = leafnode.isValueOn(pos) || leafnode.isValueOn(pos + nvo);
3100 if (isActive && (isInsideValue(lhsAcc.get(pos), iso) !=
3101 isInsideValue(lhsAcc.get((pos + nvo)), iso))) {
3102 edgeAcc.set(leafnode.offsetToGlobalCoord(pos));
3111 template<
typename LeafNodeT,
typename TreeAcc,
typename VoxelEdgeAcc>
3113 evalExternalVoxelEdges(VoxelEdgeAcc& edgeAcc,
3115 const LeafNodeT& lhsNode,
3116 const LeafNodeVoxelOffsets& voxels,
3117 const typename LeafNodeT::ValueType iso)
3119 const std::vector<Index>* lhsOffsets = &voxels.maxX();
3120 const std::vector<Index>* rhsOffsets = &voxels.minX();
3121 Coord ijk = lhsNode.origin();
3123 if (VoxelEdgeAcc::AXIS == 0) {
3124 ijk[0] += LeafNodeT::DIM;
3125 }
else if (VoxelEdgeAcc::AXIS == 1) {
3126 ijk[1] += LeafNodeT::DIM;
3127 lhsOffsets = &voxels.maxY();
3128 rhsOffsets = &voxels.minY();
3129 }
else if (VoxelEdgeAcc::AXIS == 2) {
3130 ijk[2] += LeafNodeT::DIM;
3131 lhsOffsets = &voxels.maxZ();
3132 rhsOffsets = &voxels.minZ();
3135 typename LeafNodeT::ValueType value;
3136 const LeafNodeT* rhsNodePt = acc.probeConstLeaf(ijk);
3138 const LeafBufferAccessor<LeafNodeT> lhsAcc(lhsNode);
3141 const LeafBufferAccessor<LeafNodeT> rhsAcc(*rhsNodePt);
3143 for (
size_t n = 0, N = lhsOffsets->size(); n < N; ++n) {
3144 const Index& pos = (*lhsOffsets)[n];
3145 bool isActive = lhsNode.isValueOn(pos) || rhsNodePt->isValueOn((*rhsOffsets)[n]);
3146 if (isActive && (isInsideValue(lhsAcc.get(pos), iso) !=
3147 isInsideValue(rhsAcc.get((*rhsOffsets)[n]), iso))) {
3148 edgeAcc.set(lhsNode.offsetToGlobalCoord(pos));
3151 }
else if (!acc.probeValue(ijk, value)) {
3152 const bool inside = isInsideValue(value, iso);
3153 for (
size_t n = 0, N = lhsOffsets->size(); n < N; ++n) {
3154 const Index& pos = (*lhsOffsets)[n];
3155 if (lhsNode.isValueOn(pos) && (inside != isInsideValue(lhsAcc.get(pos), iso))) {
3156 edgeAcc.set(lhsNode.offsetToGlobalCoord(pos));
3166 template<
typename LeafNodeT,
typename TreeAcc,
typename VoxelEdgeAcc>
3168 evalExternalVoxelEdgesInv(VoxelEdgeAcc& edgeAcc,
3170 const LeafNodeT& leafnode,
3171 const LeafNodeVoxelOffsets& voxels,
3172 const typename LeafNodeT::ValueType iso)
3174 Coord ijk = leafnode.origin();
3175 if (VoxelEdgeAcc::AXIS == 0) --ijk[0];
3176 else if (VoxelEdgeAcc::AXIS == 1) --ijk[1];
3177 else if (VoxelEdgeAcc::AXIS == 2) --ijk[2];
3179 typename LeafNodeT::ValueType value;
3180 if (!acc.probeConstLeaf(ijk) && !acc.probeValue(ijk, value)) {
3182 const std::vector<Index>* offsets = &voxels.internalNeighborsX();
3183 if (VoxelEdgeAcc::AXIS == 1) offsets = &voxels.internalNeighborsY();
3184 else if (VoxelEdgeAcc::AXIS == 2) offsets = &voxels.internalNeighborsZ();
3186 const LeafBufferAccessor<LeafNodeT> lhsAcc(leafnode);
3188 const bool inside = isInsideValue(value, iso);
3189 for (
size_t n = 0, N = offsets->size(); n < N; ++n) {
3191 const Index& pos = (*offsets)[n];
3192 if (leafnode.isValueOn(pos)
3193 && (inside != isInsideValue(lhsAcc.get(pos), iso)))
3195 ijk = leafnode.offsetToGlobalCoord(pos);
3196 if (VoxelEdgeAcc::AXIS == 0) --ijk[0];
3197 else if (VoxelEdgeAcc::AXIS == 1) --ijk[1];
3198 else if (VoxelEdgeAcc::AXIS == 2) --ijk[2];
3207 template<
typename InputTreeType>
3208 struct IdentifyIntersectingVoxels
3210 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
3211 using InputValueType =
typename InputLeafNodeType::ValueType;
3213 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
3215 IdentifyIntersectingVoxels(
3216 const InputTreeType& inputTree,
3217 const std::vector<const InputLeafNodeType*>& inputLeafNodes,
3218 const LeafNodeVoxelOffsets& offsets,
3219 BoolTreeType& intersectionTree,
3220 InputValueType iso);
3222 IdentifyIntersectingVoxels(IdentifyIntersectingVoxels&,
tbb::split);
3223 void operator()(
const tbb::blocked_range<size_t>&);
3224 void join(
const IdentifyIntersectingVoxels& rhs) {
3225 mIntersectionAccessor.tree().merge(rhs.mIntersectionAccessor.tree());
3230 InputLeafNodeType
const *
const *
const mInputNodes;
3232 BoolTreeType mIntersectionTree;
3235 const LeafNodeVoxelOffsets& mOffsets;
3236 const InputValueType mIsovalue;
3240 template<
typename InputTreeType>
3241 IdentifyIntersectingVoxels<InputTreeType>::IdentifyIntersectingVoxels(
3242 const InputTreeType& inputTree,
3243 const std::vector<const InputLeafNodeType*>& inputLeafNodes,
3244 const LeafNodeVoxelOffsets& offsets,
3245 BoolTreeType& intersectionTree,
3247 : mInputAccessor(inputTree)
3248 , mInputNodes(inputLeafNodes.data())
3249 , mIntersectionTree(
false)
3250 , mIntersectionAccessor(intersectionTree)
3257 template<
typename InputTreeType>
3258 IdentifyIntersectingVoxels<InputTreeType>::IdentifyIntersectingVoxels(
3260 : mInputAccessor(rhs.mInputAccessor.tree())
3261 , mInputNodes(rhs.mInputNodes)
3262 , mIntersectionTree(
false)
3263 , mIntersectionAccessor(mIntersectionTree)
3264 , mOffsets(rhs.mOffsets)
3265 , mIsovalue(rhs.mIsovalue)
3270 template<
typename InputTreeType>
3272 IdentifyIntersectingVoxels<InputTreeType>::operator()(
const tbb::blocked_range<size_t>& range)
3274 VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 0> xEdgeAcc(mIntersectionAccessor);
3275 VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 1> yEdgeAcc(mIntersectionAccessor);
3276 VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 2> zEdgeAcc(mIntersectionAccessor);
3278 for (
size_t n = range.begin(); n != range.end(); ++n) {
3280 const InputLeafNodeType& node = *mInputNodes[n];
3283 evalInternalVoxelEdges(xEdgeAcc, node, mOffsets, mIsovalue);
3285 evalInternalVoxelEdges(yEdgeAcc, node, mOffsets, mIsovalue);
3287 evalInternalVoxelEdges(zEdgeAcc, node, mOffsets, mIsovalue);
3290 evalExternalVoxelEdges(xEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3292 evalExternalVoxelEdges(yEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3294 evalExternalVoxelEdges(zEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3300 evalExternalVoxelEdgesInv(xEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3302 evalExternalVoxelEdgesInv(yEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3304 evalExternalVoxelEdgesInv(zEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3309 template<
typename InputTreeType>
3311 identifySurfaceIntersectingVoxels(
3312 typename InputTreeType::template ValueConverter<bool>::Type& intersectionTree,
3313 const InputTreeType& inputTree,
3314 typename InputTreeType::ValueType isovalue)
3316 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
3318 std::vector<const InputLeafNodeType*> inputLeafNodes;
3319 inputTree.getNodes(inputLeafNodes);
3321 LeafNodeVoxelOffsets offsets;
3322 offsets.constructOffsetList<InputLeafNodeType>();
3324 IdentifyIntersectingVoxels<InputTreeType>
op(
3325 inputTree, inputLeafNodes, offsets, intersectionTree, isovalue);
3327 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, inputLeafNodes.size()), op);
3329 maskActiveTileBorders(inputTree, isovalue, intersectionTree);
3336 template<
typename InputTreeType>
3337 struct MaskIntersectingVoxels
3339 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
3340 using InputValueType =
typename InputLeafNodeType::ValueType;
3342 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
3343 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
3345 MaskIntersectingVoxels(
3346 const InputTreeType& inputTree,
3347 const std::vector<BoolLeafNodeType*>& nodes,
3348 BoolTreeType& intersectionTree,
3349 InputValueType iso);
3351 MaskIntersectingVoxels(MaskIntersectingVoxels&,
tbb::split);
3352 void operator()(
const tbb::blocked_range<size_t>&);
3353 void join(
const MaskIntersectingVoxels& rhs) {
3354 mIntersectionAccessor.tree().merge(rhs.mIntersectionAccessor.tree());
3359 BoolLeafNodeType
const *
const *
const mNodes;
3361 BoolTreeType mIntersectionTree;
3364 const InputValueType mIsovalue;
3368 template<
typename InputTreeType>
3369 MaskIntersectingVoxels<InputTreeType>::MaskIntersectingVoxels(
3370 const InputTreeType& inputTree,
3371 const std::vector<BoolLeafNodeType*>& nodes,
3372 BoolTreeType& intersectionTree,
3374 : mInputAccessor(inputTree)
3375 , mNodes(nodes.data())
3376 , mIntersectionTree(
false)
3377 , mIntersectionAccessor(intersectionTree)
3383 template<
typename InputTreeType>
3384 MaskIntersectingVoxels<InputTreeType>::MaskIntersectingVoxels(
3386 : mInputAccessor(rhs.mInputAccessor.tree())
3387 , mNodes(rhs.mNodes)
3388 , mIntersectionTree(
false)
3389 , mIntersectionAccessor(mIntersectionTree)
3390 , mIsovalue(rhs.mIsovalue)
3395 template<
typename InputTreeType>
3397 MaskIntersectingVoxels<InputTreeType>::operator()(
const tbb::blocked_range<size_t>& range)
3399 VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 0> xEdgeAcc(mIntersectionAccessor);
3400 VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 1> yEdgeAcc(mIntersectionAccessor);
3401 VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 2> zEdgeAcc(mIntersectionAccessor);
3405 for (
size_t n = range.begin(); n != range.end(); ++n) {
3407 const BoolLeafNodeType& node = *mNodes[n];
3409 for (
auto it = node.cbeginValueOn(); it; ++it) {
3411 if (!it.getValue()) {
3413 ijk = it.getCoord();
3415 const bool inside = isInsideValue(mInputAccessor.getValue(ijk), mIsovalue);
3417 if (inside != isInsideValue(mInputAccessor.getValue(ijk.offsetBy(1, 0, 0)), mIsovalue)) {
3421 if (inside != isInsideValue(mInputAccessor.getValue(ijk.offsetBy(0, 1, 0)), mIsovalue)) {
3425 if (inside != isInsideValue(mInputAccessor.getValue(ijk.offsetBy(0, 0, 1)), mIsovalue)) {
3434 template<
typename BoolTreeType>
3435 struct MaskBorderVoxels
3437 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
3439 MaskBorderVoxels(
const BoolTreeType& maskTree,
3440 const std::vector<BoolLeafNodeType*>& maskNodes,
3441 BoolTreeType& borderTree)
3442 : mMaskTree(&maskTree)
3443 , mMaskNodes(maskNodes.data())
3444 , mTmpBorderTree(
false)
3445 , mBorderTree(&borderTree) {}
3447 MaskBorderVoxels(MaskBorderVoxels& rhs,
tbb::split)
3448 : mMaskTree(rhs.mMaskTree)
3449 , mMaskNodes(rhs.mMaskNodes)
3450 , mTmpBorderTree(
false)
3451 , mBorderTree(&mTmpBorderTree) {}
3453 void join(MaskBorderVoxels& rhs) { mBorderTree->merge(*rhs.mBorderTree); }
3455 void operator()(
const tbb::blocked_range<size_t>& range)
3461 for (
size_t n = range.begin(); n != range.end(); ++n) {
3463 const BoolLeafNodeType& node = *mMaskNodes[n];
3465 for (
auto it = node.cbeginValueOn(); it; ++it) {
3467 ijk = it.getCoord();
3469 const bool lhs = it.getValue();
3472 bool isEdgeVoxel =
false;
3475 isEdgeVoxel = (maskAcc.
probeValue(ijk, rhs) && lhs != rhs);
3478 isEdgeVoxel = isEdgeVoxel || (maskAcc.
probeValue(ijk, rhs) && lhs != rhs);
3481 isEdgeVoxel = isEdgeVoxel || (maskAcc.
probeValue(ijk, rhs) && lhs != rhs);
3484 isEdgeVoxel = isEdgeVoxel || (maskAcc.
probeValue(ijk, rhs) && lhs != rhs);
3488 isEdgeVoxel = 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);
3505 BoolTreeType
const *
const mMaskTree;
3506 BoolLeafNodeType
const *
const *
const mMaskNodes;
3508 BoolTreeType mTmpBorderTree;
3509 BoolTreeType *
const mBorderTree;
3513 template<
typename BoolTreeType>
3514 struct SyncMaskValues
3516 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
3518 SyncMaskValues(
const std::vector<BoolLeafNodeType*>& nodes,
const BoolTreeType& mask)
3519 : mNodes(nodes.data())
3520 , mMaskTree(&mask) {}
3522 void operator()(
const tbb::blocked_range<size_t>& range)
const 3524 using ValueOnIter =
typename BoolLeafNodeType::ValueOnIter;
3528 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
3530 BoolLeafNodeType& node = *mNodes[n];
3532 const BoolLeafNodeType * maskNode = maskTreeAcc.
probeConstLeaf(node.origin());
3535 for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3536 const Index pos = it.pos();
3537 if (maskNode->getValue(pos)) {
3538 node.setValueOnly(pos,
true);
3546 BoolLeafNodeType *
const *
const mNodes;
3547 BoolTreeType
const *
const mMaskTree;
3554 template<
typename BoolTreeType>
3557 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
3559 MaskSurface(
const std::vector<BoolLeafNodeType*>& nodes,
3560 const BoolTreeType& mask,
3564 : mNodes(nodes.data())
3566 , mInputTransform(inputTransform)
3567 , mMaskTransform(maskTransform)
3568 , mMatchingTransforms(mInputTransform == mMaskTransform)
3569 , mInvertMask(invert) {}
3571 void operator()(
const tbb::blocked_range<size_t>& range)
const 3573 using ValueOnIter =
typename BoolLeafNodeType::ValueOnIter;
3577 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
3579 BoolLeafNodeType& node = *mNodes[n];
3581 if (mMatchingTransforms) {
3583 const BoolLeafNodeType * maskNode = maskTreeAcc.
probeConstLeaf(node.origin());
3587 for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3588 const Index pos = it.pos();
3589 if (maskNode->isValueOn(pos) == mInvertMask) {
3590 node.setValueOnly(pos,
true);
3596 if (maskTreeAcc.
isValueOn(node.origin()) == mInvertMask) {
3597 for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3598 node.setValueOnly(it.pos(),
true);
3608 for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3610 ijk = mMaskTransform.worldToIndexCellCentered(
3611 mInputTransform.indexToWorld(it.getCoord()));
3613 if (maskTreeAcc.
isValueOn(ijk) == mInvertMask) {
3614 node.setValueOnly(it.pos(),
true);
3623 BoolLeafNodeType *
const *
const mNodes;
3624 BoolTreeType
const *
const mMaskTree;
3627 const bool mMatchingTransforms;
3628 const bool mInvertMask;
3632 template<
typename InputGr
idType>
3635 typename InputGridType::TreeType::template ValueConverter<bool>::Type& intersectionTree,
3636 typename InputGridType::TreeType::template ValueConverter<bool>::Type& borderTree,
3637 const InputGridType& inputGrid,
3639 const bool invertMask,
3640 const typename InputGridType::ValueType isovalue)
3642 using InputTreeType =
typename InputGridType::TreeType;
3643 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
3644 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
3647 if (!maskGrid)
return;
3648 if (maskGrid->type() != BoolGridType::gridType())
return;
3651 const InputTreeType& inputTree = inputGrid.tree();
3653 const BoolGridType * surfaceMask =
static_cast<const BoolGridType*
>(maskGrid.get());
3655 const BoolTreeType& maskTree = surfaceMask->tree();
3660 std::vector<BoolLeafNodeType*> intersectionLeafNodes;
3661 intersectionTree.getNodes(intersectionLeafNodes);
3663 const tbb::blocked_range<size_t> intersectionRange(0, intersectionLeafNodes.size());
3665 tbb::parallel_for(intersectionRange,
3666 MaskSurface<BoolTreeType>(
3667 intersectionLeafNodes, maskTree, transform, maskTransform, invertMask));
3672 MaskBorderVoxels<BoolTreeType> borderOp(
3673 intersectionTree, intersectionLeafNodes, borderTree);
3674 tbb::parallel_reduce(intersectionRange, borderOp);
3679 BoolTreeType tmpIntersectionTree(
false);
3681 MaskIntersectingVoxels<InputTreeType>
op(
3682 inputTree, intersectionLeafNodes, tmpIntersectionTree, isovalue);
3684 tbb::parallel_reduce(intersectionRange, op);
3686 std::vector<BoolLeafNodeType*> tmpIntersectionLeafNodes;
3687 tmpIntersectionTree.getNodes(tmpIntersectionLeafNodes);
3689 tbb::parallel_for(tbb::blocked_range<size_t>(0, tmpIntersectionLeafNodes.size()),
3690 SyncMaskValues<BoolTreeType>(tmpIntersectionLeafNodes, intersectionTree));
3692 intersectionTree.clear();
3693 intersectionTree.merge(tmpIntersectionTree);
3700 template<
typename InputTreeType>
3701 struct ComputeAuxiliaryData
3703 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
3704 using InputValueType =
typename InputLeafNodeType::ValueType;
3708 using Int16TreeType =
typename InputTreeType::template ValueConverter<Int16>::Type;
3709 using Index32TreeType =
typename InputTreeType::template ValueConverter<Index32>::Type;
3712 ComputeAuxiliaryData(
const InputTreeType& inputTree,
3713 const std::vector<const BoolLeafNodeType*>& intersectionLeafNodes,
3714 Int16TreeType& signFlagsTree,
3715 Index32TreeType& pointIndexTree,
3716 InputValueType iso);
3718 ComputeAuxiliaryData(ComputeAuxiliaryData&,
tbb::split);
3719 void operator()(
const tbb::blocked_range<size_t>&);
3720 void join(
const ComputeAuxiliaryData& rhs) {
3721 mSignFlagsAccessor.tree().merge(rhs.mSignFlagsAccessor.tree());
3722 mPointIndexAccessor.tree().merge(rhs.mPointIndexAccessor.tree());
3727 BoolLeafNodeType
const *
const *
const mIntersectionNodes;
3729 Int16TreeType mSignFlagsTree;
3731 Index32TreeType mPointIndexTree;
3734 const InputValueType mIsovalue;
3738 template<
typename InputTreeType>
3739 ComputeAuxiliaryData<InputTreeType>::ComputeAuxiliaryData(
3740 const InputTreeType& inputTree,
3741 const std::vector<const BoolLeafNodeType*>& intersectionLeafNodes,
3742 Int16TreeType& signFlagsTree,
3743 Index32TreeType& pointIndexTree,
3745 : mInputAccessor(inputTree)
3746 , mIntersectionNodes(intersectionLeafNodes.data())
3748 , mSignFlagsAccessor(signFlagsTree)
3750 , mPointIndexAccessor(pointIndexTree)
3757 template<
typename InputTreeType>
3758 ComputeAuxiliaryData<InputTreeType>::ComputeAuxiliaryData(ComputeAuxiliaryData& rhs,
tbb::split)
3759 : mInputAccessor(rhs.mInputAccessor.tree())
3760 , mIntersectionNodes(rhs.mIntersectionNodes)
3762 , mSignFlagsAccessor(mSignFlagsTree)
3764 , mPointIndexAccessor(mPointIndexTree)
3765 , mIsovalue(rhs.mIsovalue)
3770 template<
typename InputTreeType>
3772 ComputeAuxiliaryData<InputTreeType>::operator()(
const tbb::blocked_range<size_t>& range)
3774 using Int16LeafNodeType =
typename Int16TreeType::LeafNodeType;
3777 std::array<InputValueType, 8> cellVertexValues;
3778 std::unique_ptr<Int16LeafNodeType> signsNodePt(
nullptr);
3780 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
3782 const BoolLeafNodeType& maskNode = *mIntersectionNodes[n];
3783 const Coord& origin = maskNode.origin();
3785 const InputLeafNodeType* leafPt = mInputAccessor.probeConstLeaf(origin);
3787 if (!signsNodePt.get()) signsNodePt.reset(
new Int16LeafNodeType(origin, 0));
3788 else signsNodePt->setOrigin(origin);
3790 bool updatedNode =
false;
3792 for (
auto it = maskNode.cbeginValueOn(); it; ++it) {
3794 const Index pos = it.pos();
3795 ijk = BoolLeafNodeType::offsetToLocalCoord(pos);
3797 const bool inclusiveCell = leafPt && isInternalLeafCoord<InputLeafNodeType>(ijk);
3799 if (inclusiveCell) getCellVertexValues(*leafPt, pos, cellVertexValues);
3800 else getCellVertexValues(mInputAccessor, origin + ijk, cellVertexValues);
3802 uint8_t signFlags = computeSignFlags(cellVertexValues, mIsovalue);
3804 if (signFlags != 0 && signFlags != 0xFF) {
3806 const bool inside = signFlags & 0x1;
3808 int edgeFlags = inside ? INSIDE : 0;
3810 if (!it.getValue()) {
3811 edgeFlags |= inside != ((signFlags & 0x02) != 0) ? XEDGE : 0;
3812 edgeFlags |= inside != ((signFlags & 0x10) != 0) ? YEDGE : 0;
3813 edgeFlags |= inside != ((signFlags & 0x08) != 0) ? ZEDGE : 0;
3816 const uint8_t ambiguousCellFlags = sAmbiguousFace[signFlags];
3817 if (ambiguousCellFlags != 0) {
3818 correctCellSigns(signFlags, ambiguousCellFlags, mInputAccessor,
3819 origin + ijk, mIsovalue);
3822 edgeFlags |= int(signFlags);
3824 signsNodePt->setValueOn(pos,
Int16(edgeFlags));
3830 typename Index32TreeType::LeafNodeType* idxNode =
3831 mPointIndexAccessor.touchLeaf(origin);
3832 idxNode->topologyUnion(*signsNodePt);
3835 auto*
const idxData = idxNode->buffer().data();
3836 for (
auto it = idxNode->beginValueOn(); it; ++it) {
3837 idxData[it.pos()] = 0;
3840 mSignFlagsAccessor.addLeaf(signsNodePt.release());
3846 template<
typename InputTreeType>
3848 computeAuxiliaryData(
3849 typename InputTreeType::template ValueConverter<Int16>::Type& signFlagsTree,
3850 typename InputTreeType::template ValueConverter<Index32>::Type& pointIndexTree,
3851 const typename InputTreeType::template ValueConverter<bool>::Type& intersectionTree,
3852 const InputTreeType& inputTree,
3853 typename InputTreeType::ValueType isovalue)
3855 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
3856 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
3858 std::vector<const BoolLeafNodeType*> intersectionLeafNodes;
3859 intersectionTree.getNodes(intersectionLeafNodes);
3861 ComputeAuxiliaryData<InputTreeType>
op(
3862 inputTree, intersectionLeafNodes, signFlagsTree, pointIndexTree, isovalue);
3864 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, intersectionLeafNodes.size()), op);
3871 template<Index32 LeafNodeLog2Dim>
3872 struct LeafNodePointCount
3876 LeafNodePointCount(
const std::vector<Int16LeafNodeType*>& leafNodes,
3877 std::unique_ptr<
Index32[]>& leafNodeCount)
3878 : mLeafNodes(leafNodes.data())
3879 , mData(leafNodeCount.get()) {}
3881 void operator()(
const tbb::blocked_range<size_t>& range)
const {
3883 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
3887 Int16 const * p = mLeafNodes[n]->buffer().data();
3888 Int16 const *
const endP = p + Int16LeafNodeType::SIZE;
3891 count +=
Index32(sEdgeGroupTable[(SIGNS & *p)][0]);
3900 Int16LeafNodeType *
const *
const mLeafNodes;
3905 template<
typename Po
intIndexLeafNode>
3906 struct AdaptiveLeafNodePointCount
3910 AdaptiveLeafNodePointCount(
const std::vector<PointIndexLeafNode*>& pointIndexNodes,
3911 const std::vector<Int16LeafNodeType*>& signDataNodes,
3912 std::unique_ptr<
Index32[]>& leafNodeCount)
3913 : mPointIndexNodes(pointIndexNodes.data())
3914 , mSignDataNodes(signDataNodes.data())
3915 , mData(leafNodeCount.get()) {}
3917 void operator()(
const tbb::blocked_range<size_t>& range)
const 3921 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
3924 const Int16LeafNodeType& signNode = *mSignDataNodes[n];
3928 std::set<IndexType> uniqueRegions;
3935 count += size_t(sEdgeGroupTable[(SIGNS & signNode.getValue(it.pos()))][0]);
3937 uniqueRegions.insert(
id);
3941 mData[n] =
Index32(count + uniqueRegions.size());
3947 Int16LeafNodeType
const *
const *
const mSignDataNodes;
3952 template<
typename Po
intIndexLeafNode>
3957 MapPoints(
const std::vector<PointIndexLeafNode*>& pointIndexNodes,
3958 const std::vector<Int16LeafNodeType*>& signDataNodes,
3959 std::unique_ptr<
Index32[]>& leafNodeCount)
3960 : mPointIndexNodes(pointIndexNodes.data())
3961 , mSignDataNodes(signDataNodes.data())
3962 , mData(leafNodeCount.get()) {}
3964 void operator()(
const tbb::blocked_range<size_t>& range)
const {
3966 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
3968 const Int16LeafNodeType& signNode = *mSignDataNodes[n];
3971 Index32 pointOffset = mData[n];
3974 const Index pos = it.pos();
3976 const int signs = SIGNS & int(signNode.getValue(pos));
3977 pointOffset +=
Index32(sEdgeGroupTable[signs][0]);
3984 Int16LeafNodeType
const *
const *
const mSignDataNodes;
3989 template<
typename TreeType,
typename PrimBuilder>
3990 struct ComputePolygons
3992 using Int16TreeType =
typename TreeType::template ValueConverter<Int16>::Type;
3993 using Int16LeafNodeType =
typename Int16TreeType::LeafNodeType;
3995 using Index32TreeType =
typename TreeType::template ValueConverter<Index32>::Type;
3996 using Index32LeafNodeType =
typename Index32TreeType::LeafNodeType;
3999 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
4000 const Int16TreeType& signFlagsTree,
4001 const Index32TreeType& idxTree,
4003 bool invertSurfaceOrientation);
4005 void setRefSignTree(
const Int16TreeType * r) { mRefSignFlagsTree = r; }
4007 void operator()(
const tbb::blocked_range<size_t>&)
const;
4010 Int16LeafNodeType *
const *
const mSignFlagsLeafNodes;
4011 Int16TreeType
const *
const mSignFlagsTree;
4012 Int16TreeType
const * mRefSignFlagsTree;
4013 Index32TreeType
const *
const mIndexTree;
4015 bool const mInvertSurfaceOrientation;
4019 template<
typename TreeType,
typename PrimBuilder>
4020 ComputePolygons<TreeType, PrimBuilder>::ComputePolygons(
4021 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
4022 const Int16TreeType& signFlagsTree,
4023 const Index32TreeType& idxTree,
4025 bool invertSurfaceOrientation)
4026 : mSignFlagsLeafNodes(signFlagsLeafNodes.data())
4027 , mSignFlagsTree(&signFlagsTree)
4028 , mRefSignFlagsTree(
nullptr)
4029 , mIndexTree(&idxTree)
4030 , mPolygonPoolList(&polygons)
4031 , mInvertSurfaceOrientation(invertSurfaceOrientation)
4035 template<
typename InputTreeType,
typename PrimBuilder>
4037 ComputePolygons<InputTreeType, PrimBuilder>::operator()(
const tbb::blocked_range<size_t>& range)
const 4040 Int16ValueAccessor signAcc(*mSignFlagsTree);
4044 const bool invertSurfaceOrientation = mInvertSurfaceOrientation;
4051 std::unique_ptr<Int16ValueAccessor> refSignAcc;
4052 if (mRefSignFlagsTree) refSignAcc.reset(
new Int16ValueAccessor(*mRefSignFlagsTree));
4054 for (
size_t n = range.begin(); n != range.end(); ++n) {
4056 const Int16LeafNodeType& node = *mSignFlagsLeafNodes[n];
4057 origin = node.origin();
4061 typename Int16LeafNodeType::ValueOnCIter iter = node.cbeginValueOn();
4062 for (; iter; ++iter) {
4063 if (iter.getValue() & XEDGE) ++edgeCount;
4064 if (iter.getValue() & YEDGE) ++edgeCount;
4065 if (iter.getValue() & ZEDGE) ++edgeCount;
4068 if (edgeCount == 0)
continue;
4070 mesher.init(edgeCount, (*mPolygonPoolList)[n]);
4072 const Int16LeafNodeType *signleafPt = signAcc.probeConstLeaf(origin);
4073 const Index32LeafNodeType *idxLeafPt = idxAcc.
probeConstLeaf(origin);
4075 if (!signleafPt || !idxLeafPt)
continue;
4077 const Int16LeafNodeType *refSignLeafPt =
nullptr;
4078 if (refSignAcc) refSignLeafPt = refSignAcc->probeConstLeaf(origin);
4082 for (iter = node.cbeginValueOn(); iter; ++iter) {
4083 ijk = iter.getCoord();
4085 const Int16 flags = iter.getValue();
4086 if (!(flags & 0xE00))
continue;
4089 if (refSignLeafPt) {
4090 refFlags = refSignLeafPt->getValue(iter.pos());
4093 const uint8_t cell = uint8_t(SIGNS & flags);
4095 if (sEdgeGroupTable[cell][0] > 1) {
4096 offsets[0] = (sEdgeGroupTable[cell][1] - 1);
4097 offsets[1] = (sEdgeGroupTable[cell][9] - 1);
4098 offsets[2] = (sEdgeGroupTable[cell][4] - 1);
4104 if (ijk[0] > origin[0] && ijk[1] > origin[1] && ijk[2] > origin[2]) {
4105 constructPolygons(invertSurfaceOrientation,
4106 flags, refFlags, offsets, ijk, *signleafPt, *idxLeafPt, mesher);
4108 constructPolygons(invertSurfaceOrientation,
4109 flags, refFlags, offsets, ijk, signAcc, idxAcc, mesher);
4122 template<
typename T>
4125 CopyArray(T * outputArray,
const T * inputArray,
size_t outputOffset = 0)
4126 : mOutputArray(outputArray), mInputArray(inputArray), mOutputOffset(outputOffset)
4130 void operator()(
const tbb::blocked_range<size_t>& inputArrayRange)
const 4132 const size_t offset = mOutputOffset;
4133 for (
size_t n = inputArrayRange.begin(), N = inputArrayRange.end(); n < N; ++n) {
4134 mOutputArray[offset + n] = mInputArray[n];
4139 T *
const mOutputArray;
4140 T
const *
const mInputArray;
4141 size_t const mOutputOffset;
4145 struct FlagAndCountQuadsToSubdivide
4148 const std::vector<uint8_t>& pointFlags,
4150 std::unique_ptr<
unsigned[]>& numQuadsToDivide)
4151 : mPolygonPoolList(&polygons)
4152 , mPointFlags(pointFlags.data())
4153 , mPoints(points.get())
4154 , mNumQuadsToDivide(numQuadsToDivide.get())
4158 void operator()(
const tbb::blocked_range<size_t>& range)
const 4160 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
4167 for (
size_t i = 0, I = polygons.
numQuads(); i < I; ++i) {
4175 const bool edgePoly = mPointFlags[quad[0]] || mPointFlags[quad[1]]
4176 || mPointFlags[quad[2]] || mPointFlags[quad[3]];
4178 if (!edgePoly)
continue;
4180 const Vec3s& p0 = mPoints[quad[0]];
4181 const Vec3s& p1 = mPoints[quad[1]];
4182 const Vec3s& p2 = mPoints[quad[2]];
4183 const Vec3s& p3 = mPoints[quad[3]];
4185 if (!isPlanarQuad(p0, p1, p2, p3, 1e-6f)) {
4192 mNumQuadsToDivide[n] = count;
4198 uint8_t
const *
const mPointFlags;
4199 Vec3s const *
const mPoints;
4200 unsigned *
const mNumQuadsToDivide;
4204 struct SubdivideQuads
4210 std::unique_ptr<
unsigned[]>& numQuadsToDivide,
4211 std::unique_ptr<
unsigned[]>& centroidOffsets)
4212 : mPolygonPoolList(&polygons)
4213 , mPoints(points.get())
4214 , mCentroids(centroids.get())
4215 , mNumQuadsToDivide(numQuadsToDivide.get())
4216 , mCentroidOffsets(centroidOffsets.get())
4217 , mPointCount(pointCount)
4221 void operator()(
const tbb::blocked_range<size_t>& range)
const 4223 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
4227 const size_t nonplanarCount = size_t(mNumQuadsToDivide[n]);
4229 if (nonplanarCount > 0) {
4235 size_t offset = mCentroidOffsets[n];
4237 size_t triangleIdx = 0;
4239 for (
size_t i = 0, I = polygons.
numQuads(); i < I; ++i) {
4241 const char quadFlags = polygons.
quadFlags(i);
4244 unsigned newPointIdx = unsigned(offset + mPointCount);
4248 mCentroids[offset] = (mPoints[quad[0]] + mPoints[quad[1]] +
4249 mPoints[quad[2]] + mPoints[quad[3]]) * 0.25f;
4256 triangle[0] = quad[0];
4257 triangle[1] = newPointIdx;
4258 triangle[2] = quad[3];
4268 triangle[0] = quad[0];
4269 triangle[1] = quad[1];
4270 triangle[2] = newPointIdx;
4280 triangle[0] = quad[1];
4281 triangle[1] = quad[2];
4282 triangle[2] = newPointIdx;
4293 triangle[0] = quad[2];
4294 triangle[1] = quad[3];
4295 triangle[2] = newPointIdx;
4306 for (
size_t i = 0, I = polygons.
numTriangles(); i < I; ++i) {
4313 for (
size_t i = 0, I = polygons.
numQuads(); i < I; ++i) {
4317 tmpPolygons.
quad(quadIdx) = quad;
4323 polygons.
copy(tmpPolygons);
4330 Vec3s const *
const mPoints;
4331 Vec3s *
const mCentroids;
4332 unsigned *
const mNumQuadsToDivide;
4333 unsigned *
const mCentroidOffsets;
4334 size_t const mPointCount;
4339 subdivideNonplanarSeamLineQuads(
4341 size_t polygonPoolListSize,
4343 size_t& pointListSize,
4344 std::vector<uint8_t>& pointFlags)
4346 const tbb::blocked_range<size_t> polygonPoolListRange(0, polygonPoolListSize);
4348 std::unique_ptr<unsigned[]> numQuadsToDivide(
new unsigned[polygonPoolListSize]);
4350 tbb::parallel_for(polygonPoolListRange,
4351 FlagAndCountQuadsToSubdivide(polygonPoolList, pointFlags, pointList, numQuadsToDivide));
4353 std::unique_ptr<unsigned[]> centroidOffsets(
new unsigned[polygonPoolListSize]);
4355 size_t centroidCount = 0;
4359 for (
size_t n = 0, N = polygonPoolListSize; n < N; ++n) {
4360 centroidOffsets[n] = sum;
4361 sum += numQuadsToDivide[n];
4363 centroidCount = size_t(sum);
4366 std::unique_ptr<Vec3s[]> centroidList(
new Vec3s[centroidCount]);
4368 tbb::parallel_for(polygonPoolListRange,
4369 SubdivideQuads(polygonPoolList, pointList, pointListSize,
4370 centroidList, numQuadsToDivide, centroidOffsets));
4372 if (centroidCount > 0) {
4374 const size_t newPointListSize = centroidCount + pointListSize;
4376 std::unique_ptr<openvdb::Vec3s[]> newPointList(
new openvdb::Vec3s[newPointListSize]);
4378 tbb::parallel_for(tbb::blocked_range<size_t>(0, pointListSize),
4379 CopyArray<Vec3s>(newPointList.get(), pointList.get()));
4381 tbb::parallel_for(tbb::blocked_range<size_t>(0, newPointListSize - pointListSize),
4382 CopyArray<Vec3s>(newPointList.get(), centroidList.get(), pointListSize));
4384 pointListSize = newPointListSize;
4385 pointList.swap(newPointList);
4386 pointFlags.resize(pointListSize, 0);
4391 struct ReviseSeamLineFlags
4394 const std::vector<uint8_t>& pointFlags)
4395 : mPolygonPoolList(&polygons)
4396 , mPointFlags(pointFlags.data())
4400 void operator()(
const tbb::blocked_range<size_t>& range)
const 4402 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
4406 for (
size_t i = 0, I = polygons.
numQuads(); i < I; ++i) {
4414 const bool hasSeamLinePoint =
4415 mPointFlags[verts[0]] || mPointFlags[verts[1]] ||
4416 mPointFlags[verts[2]] || mPointFlags[verts[3]];
4418 if (!hasSeamLinePoint) {
4419 flags &= ~POLYFLAG_FRACTURE_SEAM;
4424 for (
size_t i = 0, I = polygons.
numTriangles(); i < I; ++i) {
4432 const bool hasSeamLinePoint =
4433 mPointFlags[verts[0]] || mPointFlags[verts[1]] || mPointFlags[verts[2]];
4435 if (!hasSeamLinePoint) {
4436 flags &= ~POLYFLAG_FRACTURE_SEAM;
4447 uint8_t
const *
const mPointFlags;
4452 reviseSeamLineFlags(
PolygonPoolList& polygonPoolList,
size_t polygonPoolListSize,
4453 std::vector<uint8_t>& pointFlags)
4455 tbb::parallel_for(tbb::blocked_range<size_t>(0, polygonPoolListSize),
4456 ReviseSeamLineFlags(polygonPoolList, pointFlags));
4463 template<
typename InputTreeType>
4464 struct MaskDisorientedTrianglePoints
4466 MaskDisorientedTrianglePoints(
const InputTreeType& inputTree,
const PolygonPoolList& polygons,
4467 const PointList& pointList, std::unique_ptr<uint8_t[]>& pointMask,
4469 : mInputTree(&inputTree)
4470 , mPolygonPoolList(&polygons)
4471 , mPointList(&pointList)
4472 , mPointMask(pointMask.get())
4473 , mTransform(transform)
4474 , mInvertSurfaceOrientation(invertSurfaceOrientation)
4478 void operator()(
const tbb::blocked_range<size_t>& range)
const 4480 using ValueType =
typename InputTreeType::LeafNodeType::ValueType;
4483 Vec3s centroid, normal;
4486 const bool invertGradientDir = mInvertSurfaceOrientation || isBoolValue<ValueType>();
4488 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
4490 const PolygonPool& polygons = (*mPolygonPoolList)[n];
4492 for (
size_t i = 0, I = polygons.
numTriangles(); i < I; ++i) {
4496 const Vec3s& v0 = (*mPointList)[verts[0]];
4497 const Vec3s& v1 = (*mPointList)[verts[1]];
4498 const Vec3s& v2 = (*mPointList)[verts[2]];
4500 normal = (v2 - v0).cross((v1 - v0));
4503 centroid = (v0 + v1 + v2) * (1.0f / 3.0f);
4504 ijk = mTransform.worldToIndexCellCentered(centroid);
4509 if (invertGradientDir) {
4514 if (dir.dot(normal) < -0.5f) {
4519 mPointMask[verts[0]] = 1;
4520 mPointMask[verts[1]] = 1;
4521 mPointMask[verts[2]] = 1;
4530 InputTreeType
const *
const mInputTree;
4533 uint8_t *
const mPointMask;
4535 bool const mInvertSurfaceOrientation;
4539 template<
typename InputTree>
4541 relaxDisorientedTriangles(
4542 bool invertSurfaceOrientation,
4543 const InputTree& inputTree,
4546 size_t polygonPoolListSize,
4548 const size_t pointListSize)
4550 const tbb::blocked_range<size_t> polygonPoolListRange(0, polygonPoolListSize);
4552 std::unique_ptr<uint8_t[]> pointMask(
new uint8_t[pointListSize]);
4553 fillArray(pointMask.get(), uint8_t(0), pointListSize);
4555 tbb::parallel_for(polygonPoolListRange,
4556 MaskDisorientedTrianglePoints<InputTree>(
4557 inputTree, polygonPoolList, pointList, pointMask, transform, invertSurfaceOrientation));
4559 std::unique_ptr<uint8_t[]> pointUpdates(
new uint8_t[pointListSize]);
4560 fillArray(pointUpdates.get(), uint8_t(0), pointListSize);
4562 std::unique_ptr<Vec3s[]> newPoints(
new Vec3s[pointListSize]);
4563 fillArray(newPoints.get(),
Vec3s(0.0f, 0.0f, 0.0f), pointListSize);
4565 for (
size_t n = 0, N = polygonPoolListSize; n < N; ++n) {
4569 for (
size_t i = 0, I = polygons.
numQuads(); i < I; ++i) {
4572 for (
int v = 0; v < 4; ++v) {
4574 const unsigned pointIdx = verts[v];
4576 if (pointMask[pointIdx] == 1) {
4578 newPoints[pointIdx] +=
4579 pointList[verts[0]] + pointList[verts[1]] +
4580 pointList[verts[2]] + pointList[verts[3]];
4582 pointUpdates[pointIdx] = uint8_t(pointUpdates[pointIdx] + 4);
4587 for (
size_t i = 0, I = polygons.
numTriangles(); i < I; ++i) {
4590 for (
int v = 0; v < 3; ++v) {
4592 const unsigned pointIdx = verts[v];
4594 if (pointMask[pointIdx] == 1) {
4595 newPoints[pointIdx] +=
4596 pointList[verts[0]] + pointList[verts[1]] + pointList[verts[2]];
4598 pointUpdates[pointIdx] = uint8_t(pointUpdates[pointIdx] + 3);
4604 for (
size_t n = 0, N = pointListSize; n < N; ++n) {
4605 if (pointUpdates[n] > 0) {
4606 const double weight = 1.0 / double(pointUpdates[n]);
4607 pointList[n] = newPoints[n] * float(weight);
4613 template<
typename Gr
idType>
4617 std::vector<Vec3s>& points,
4618 std::vector<Vec3I>& triangles,
4619 std::vector<Vec4I>& quads,
4622 bool relaxDisorientedTriangles)
4624 static_assert(std::is_scalar<typename GridType::ValueType>::value,
4625 "volume to mesh conversion is supported only for scalar grids");
4627 VolumeToMesh mesher(isovalue, adaptivity, relaxDisorientedTriangles);
4635 volume_to_mesh_internal::PointListCopy ptnCpy(mesher.
pointList(), points);
4636 tbb::parallel_for(tbb::blocked_range<size_t>(0, points.size()), ptnCpy);
4643 size_t numQuads = 0, numTriangles = 0;
4645 openvdb::tools::PolygonPool& polygons = polygonPoolList[n];
4646 numTriangles += polygons.numTriangles();
4647 numQuads += polygons.numQuads();
4651 triangles.resize(numTriangles);
4653 quads.resize(numQuads);
4657 size_t qIdx = 0, tIdx = 0;
4659 openvdb::tools::PolygonPool& polygons = polygonPoolList[n];
4661 for (
size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4662 quads[qIdx++] = polygons.quad(i);
4665 for (
size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
4666 triangles[tIdx++] = polygons.triangle(i);
4680 PolygonPool::PolygonPool()
4684 , mTriangles(nullptr)
4685 , mQuadFlags(nullptr)
4686 , mTriangleFlags(nullptr)
4693 : mNumQuads(numQuads)
4694 , mNumTriangles(numTriangles)
4697 , mQuadFlags(new char[mNumQuads])
4698 , mTriangleFlags(new char[mNumTriangles])
4709 for (
size_t i = 0; i < mNumQuads; ++i) {
4710 mQuads[i] = rhs.mQuads[i];
4711 mQuadFlags[i] = rhs.mQuadFlags[i];
4714 for (
size_t i = 0; i < mNumTriangles; ++i) {
4715 mTriangles[i] = rhs.mTriangles[i];
4716 mTriangleFlags[i] = rhs.mTriangleFlags[i];
4726 mQuadFlags.reset(
new char[mNumQuads]);
4734 mQuads.reset(
nullptr);
4735 mQuadFlags.reset(
nullptr);
4742 mNumTriangles = size;
4744 mTriangleFlags.reset(
new char[mNumTriangles]);
4752 mTriangles.reset(
nullptr);
4753 mTriangleFlags.reset(
nullptr);
4760 if (!(n < mNumQuads))
return false;
4765 mQuads.reset(
nullptr);
4769 std::unique_ptr<char[]> flags(
new char[n]);
4771 for (
size_t i = 0; i < n; ++i) {
4772 quads[i] = mQuads[i];
4773 flags[i] = mQuadFlags[i];
4777 mQuadFlags.swap(flags);
4789 if (!(n < mNumTriangles))
return false;
4794 mTriangles.reset(
nullptr);
4797 std::unique_ptr<openvdb::Vec3I[]> triangles(
new openvdb::Vec3I[n]);
4798 std::unique_ptr<char[]> flags(
new char[n]);
4800 for (
size_t i = 0; i < n; ++i) {
4801 triangles[i] = mTriangles[i];
4802 flags[i] = mTriangleFlags[i];
4805 mTriangles.swap(triangles);
4806 mTriangleFlags.swap(flags);
4823 , mSeamPointListSize(0)
4824 , mPolygonPoolListSize(0)
4825 , mIsovalue(isovalue)
4826 , mPrimAdaptivity(adaptivity)
4827 , mSecAdaptivity(0.0)
4829 , mSurfaceMaskGrid(
GridBase::ConstPtr())
4830 , mAdaptivityGrid(
GridBase::ConstPtr())
4831 , mAdaptivityMaskTree(
TreeBase::ConstPtr())
4834 , mInvertSurfaceMask(false)
4835 , mRelaxDisorientedTriangles(relaxDisorientedTriangles)
4836 , mQuantizedSeamPoints(nullptr)
4846 mSecAdaptivity = secAdaptivity;
4851 mSeamPointListSize = 0;
4852 mQuantizedSeamPoints.reset(
nullptr);
4859 mSurfaceMaskGrid = mask;
4860 mInvertSurfaceMask = invertMask;
4867 mAdaptivityGrid = grid;
4874 mAdaptivityMaskTree = tree;
4878 template<
typename InputGr
idType>
4884 using InputTreeType =
typename InputGridType::TreeType;
4885 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
4886 using InputValueType =
typename InputLeafNodeType::ValueType;
4890 using FloatTreeType =
typename InputTreeType::template ValueConverter<float>::Type;
4892 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
4893 using Int16TreeType =
typename InputTreeType::template ValueConverter<Int16>::Type;
4894 using Int16LeafNodeType =
typename Int16TreeType::LeafNodeType;
4895 using Index32TreeType =
typename InputTreeType::template ValueConverter<Index32>::Type;
4896 using Index32LeafNodeType =
typename Index32TreeType::LeafNodeType;
4901 mPolygonPoolListSize = 0;
4903 mPointFlags.clear();
4908 const InputValueType isovalue = InputValueType(mIsovalue);
4909 const float adaptivityThreshold = float(mPrimAdaptivity);
4910 const bool adaptive = mPrimAdaptivity > 1e-7 || mSecAdaptivity > 1e-7;
4916 const bool invertSurfaceOrientation = (!volume_to_mesh_internal::isBoolValue<InputValueType>()
4921 const InputTreeType& inputTree = inputGrid.tree();
4923 BoolTreeType intersectionTree(
false), adaptivityMask(
false);
4925 if (mAdaptivityMaskTree && mAdaptivityMaskTree->type() == BoolTreeType::treeType()) {
4926 const BoolTreeType *refAdaptivityMask=
4927 static_cast<const BoolTreeType*
>(mAdaptivityMaskTree.get());
4928 adaptivityMask.topologyUnion(*refAdaptivityMask);
4931 Int16TreeType signFlagsTree(0);
4937 volume_to_mesh_internal::identifySurfaceIntersectingVoxels(
4938 intersectionTree, inputTree, isovalue);
4940 volume_to_mesh_internal::applySurfaceMask(intersectionTree, adaptivityMask,
4941 inputGrid, mSurfaceMaskGrid, mInvertSurfaceMask, isovalue);
4943 if (intersectionTree.empty())
return;
4945 volume_to_mesh_internal::computeAuxiliaryData(
4946 signFlagsTree, pointIndexTree, intersectionTree, inputTree, isovalue);
4948 intersectionTree.clear();
4950 std::vector<Index32LeafNodeType*> pointIndexLeafNodes;
4951 pointIndexTree.getNodes(pointIndexLeafNodes);
4953 std::vector<Int16LeafNodeType*> signFlagsLeafNodes;
4954 signFlagsTree.getNodes(signFlagsLeafNodes);
4956 const tbb::blocked_range<size_t> auxiliaryLeafNodeRange(0, signFlagsLeafNodes.size());
4961 Int16TreeType* refSignFlagsTree =
nullptr;
4962 Index32TreeType* refPointIndexTree =
nullptr;
4963 InputTreeType
const* refInputTree =
nullptr;
4965 if (mRefGrid && mRefGrid->type() == InputGridType::gridType()) {
4967 const InputGridType* refGrid =
static_cast<const InputGridType*
>(mRefGrid.get());
4968 refInputTree = &refGrid->tree();
4970 if (!mRefSignTree && !mRefIdxTree) {
4974 typename Int16TreeType::Ptr refSignFlagsTreePt(
new Int16TreeType(0));
4975 typename Index32TreeType::Ptr refPointIndexTreePt(
4978 BoolTreeType refIntersectionTree(
false);
4980 volume_to_mesh_internal::identifySurfaceIntersectingVoxels(
4981 refIntersectionTree, *refInputTree, isovalue);
4983 volume_to_mesh_internal::computeAuxiliaryData(*refSignFlagsTreePt,
4984 *refPointIndexTreePt, refIntersectionTree, *refInputTree, isovalue);
4986 mRefSignTree = refSignFlagsTreePt;
4987 mRefIdxTree = refPointIndexTreePt;
4990 if (mRefSignTree && mRefIdxTree) {
4994 refSignFlagsTree =
static_cast<Int16TreeType*
>(mRefSignTree.get());
4995 refPointIndexTree =
static_cast<Index32TreeType*
>(mRefIdxTree.get());
4999 if (refSignFlagsTree && refPointIndexTree) {
5003 volume_to_mesh_internal::markSeamLineData(signFlagsTree, *refSignFlagsTree);
5005 if (mSeamPointListSize == 0) {
5009 std::vector<Int16LeafNodeType*> refSignFlagsLeafNodes;
5010 refSignFlagsTree->getNodes(refSignFlagsLeafNodes);
5012 std::unique_ptr<Index32[]> leafNodeOffsets(
5013 new Index32[refSignFlagsLeafNodes.size()]);
5015 tbb::parallel_for(tbb::blocked_range<size_t>(0, refSignFlagsLeafNodes.size()),
5016 volume_to_mesh_internal::LeafNodePointCount<Int16LeafNodeType::LOG2DIM>(
5017 refSignFlagsLeafNodes, leafNodeOffsets));
5021 for (
size_t n = 0, N = refSignFlagsLeafNodes.size(); n < N; ++n) {
5022 const Index32 tmp = leafNodeOffsets[n];
5023 leafNodeOffsets[n] = count;
5026 mSeamPointListSize = size_t(count);
5029 if (mSeamPointListSize != 0) {
5031 mQuantizedSeamPoints.reset(
new uint32_t[mSeamPointListSize]);
5033 std::memset(mQuantizedSeamPoints.get(), 0,
sizeof(uint32_t) * mSeamPointListSize);
5035 std::vector<Index32LeafNodeType*> refPointIndexLeafNodes;
5036 refPointIndexTree->getNodes(refPointIndexLeafNodes);
5038 tbb::parallel_for(tbb::blocked_range<size_t>(0, refPointIndexLeafNodes.size()),
5039 volume_to_mesh_internal::MapPoints<Index32LeafNodeType>(
5040 refPointIndexLeafNodes, refSignFlagsLeafNodes, leafNodeOffsets));
5044 if (mSeamPointListSize != 0) {
5046 tbb::parallel_for(auxiliaryLeafNodeRange,
5047 volume_to_mesh_internal::SeamLineWeights<InputTreeType>(
5048 signFlagsLeafNodes, inputTree, *refPointIndexTree, *refSignFlagsTree,
5049 mQuantizedSeamPoints.get(), isovalue));
5054 const bool referenceMeshing = refSignFlagsTree && refPointIndexTree && refInputTree;
5059 std::unique_ptr<Index32[]> leafNodeOffsets(
new Index32[signFlagsLeafNodes.size()]);
5062 volume_to_mesh_internal::MergeVoxelRegions<InputGridType> mergeOp(
5063 inputGrid, pointIndexTree, pointIndexLeafNodes, signFlagsLeafNodes,
5064 isovalue, adaptivityThreshold, invertSurfaceOrientation);
5066 if (mAdaptivityGrid && mAdaptivityGrid->type() == FloatGridType::gridType()) {
5067 const FloatGridType* adaptivityGrid =
5068 static_cast<const FloatGridType*
>(mAdaptivityGrid.get());
5069 mergeOp.setSpatialAdaptivity(*adaptivityGrid);
5072 if (!adaptivityMask.empty()) {
5073 mergeOp.setAdaptivityMask(adaptivityMask);
5076 if (referenceMeshing) {
5077 mergeOp.setRefSignFlagsData(*refSignFlagsTree,
float(mSecAdaptivity));
5080 tbb::parallel_for(auxiliaryLeafNodeRange, mergeOp);
5082 volume_to_mesh_internal::AdaptiveLeafNodePointCount<Index32LeafNodeType>
5083 op(pointIndexLeafNodes, signFlagsLeafNodes, leafNodeOffsets);
5085 tbb::parallel_for(auxiliaryLeafNodeRange, op);
5089 volume_to_mesh_internal::LeafNodePointCount<Int16LeafNodeType::LOG2DIM>
5090 op(signFlagsLeafNodes, leafNodeOffsets);
5092 tbb::parallel_for(auxiliaryLeafNodeRange, op);
5098 for (
size_t n = 0, N = signFlagsLeafNodes.size(); n < N; ++n) {
5099 const Index32 tmp = leafNodeOffsets[n];
5104 mPointListSize = size_t(pointCount);
5106 mPointFlags.clear();
5113 volume_to_mesh_internal::ComputePoints<InputTreeType>
5114 op(mPoints.get(), inputTree, pointIndexLeafNodes,
5115 signFlagsLeafNodes, leafNodeOffsets, transform, mIsovalue);
5117 if (referenceMeshing) {
5118 mPointFlags.resize(mPointListSize);
5119 op.setRefData(*refInputTree, *refPointIndexTree, *refSignFlagsTree,
5120 mQuantizedSeamPoints.get(), mPointFlags.data());
5123 tbb::parallel_for(auxiliaryLeafNodeRange,
op);
5129 mPolygonPoolListSize = signFlagsLeafNodes.size();
5130 mPolygons.reset(
new PolygonPool[mPolygonPoolListSize]);
5134 using PrimBuilder = volume_to_mesh_internal::AdaptivePrimBuilder;
5136 volume_to_mesh_internal::ComputePolygons<Int16TreeType, PrimBuilder>
5137 op(signFlagsLeafNodes, signFlagsTree, pointIndexTree,
5138 mPolygons, invertSurfaceOrientation);
5140 if (referenceMeshing) {
5141 op.setRefSignTree(refSignFlagsTree);
5144 tbb::parallel_for(auxiliaryLeafNodeRange, op);
5148 using PrimBuilder = volume_to_mesh_internal::UniformPrimBuilder;
5150 volume_to_mesh_internal::ComputePolygons<Int16TreeType, PrimBuilder>
5151 op(signFlagsLeafNodes, signFlagsTree, pointIndexTree,
5152 mPolygons, invertSurfaceOrientation);
5154 if (referenceMeshing) {
5155 op.setRefSignTree(refSignFlagsTree);
5158 tbb::parallel_for(auxiliaryLeafNodeRange, op);
5162 signFlagsTree.clear();
5163 pointIndexTree.clear();
5166 if (adaptive && mRelaxDisorientedTriangles) {
5167 volume_to_mesh_internal::relaxDisorientedTriangles(invertSurfaceOrientation,
5168 inputTree, transform, mPolygons, mPolygonPoolListSize, mPoints, mPointListSize);
5172 if (referenceMeshing) {
5173 volume_to_mesh_internal::subdivideNonplanarSeamLineQuads(
5174 mPolygons, mPolygonPoolListSize, mPoints, mPointListSize, mPointFlags);
5176 volume_to_mesh_internal::reviseSeamLineFlags(mPolygons, mPolygonPoolListSize, mPointFlags);
5185 template<
typename Gr
idType>
5188 std::vector<Vec3s>& points,
5189 std::vector<Vec3I>& triangles,
5190 std::vector<Vec4I>& quads,
5193 bool relaxDisorientedTriangles)
5195 volume_to_mesh_internal::doVolumeToMesh(grid, points, triangles, quads,
5196 isovalue, adaptivity, relaxDisorientedTriangles);
5199 template<
typename Gr
idType>
5202 std::vector<Vec3s>& points,
5203 std::vector<Vec4I>& quads,
5206 std::vector<Vec3I> triangles;
5207 volumeToMesh(grid, points, triangles, quads, isovalue, 0.0,
true);
5216 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION 5218 #ifdef OPENVDB_INSTANTIATE_VOLUMETOMESH 5222 #define _FUNCTION(TreeT) \ 5223 void volumeToMesh(const Grid<TreeT>&, std::vector<Vec3s>&, std::vector<Vec4I>&, double) 5227 #define _FUNCTION(TreeT) \ 5228 void volumeToMesh(const Grid<TreeT>&, std::vector<Vec3s>&, std::vector<Vec3I>&, std::vector<Vec4I>&, double, double, bool) 5232 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION 5239 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
Templated block class to hold specific data types and a fixed number of values determined by Log2Dim...
Definition: LeafNode.h:38
Mat3< double > Mat3d
Definition: Mat3.h:834
The Value Accessor Implementation and API methods. The majoirty of the API matches the API of a compa...
Definition: ValueAccessor.h:68
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
int getValueDepth(const Coord &xyz) const
Return the tree depth (0 = root) at which the value of voxel (x, y, z) resides, or -1 if (x...
Definition: ValueAccessor.h:516
#define OPENVDB_NUMERIC_TREE_INSTANTIATE(Function)
Definition: version.h.in:163
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: LeafNode.h:1075
int16_t Int16
Definition: Types.h:55
Mat3 transpose() const
returns transpose of this
Definition: Mat3.h:454
T & z()
Definition: Vec3.h:88
Index32 Index
Definition: Types.h:54
constexpr Index32 INVALID_IDX
Definition: Util.h:19
Definition: Exceptions.h:65
bool normalize(T eps=T(1.0e-7))
this = normalized this
Definition: Vec3.h:363
void split(ContainerT &out, const std::string &in, const char delim)
Definition: Name.h:43
Vec3< float > Vec3s
Definition: Vec3.h:664
bool probeValue(const Coord &xyz, ValueType &value) const
Return the active state of the value at a given coordinate as well as its value.
Definition: ValueAccessor.h:492
#define OPENVDB_ASSERT(X)
Definition: Assert.h:41
Base class for typed trees.
Definition: Tree.h:37
Index64 pointCount(const PointDataTreeT &tree, const FilterT &filter=NullFilter(), const bool inCoreOnly=false, const bool threaded=true)
Count the total number of points in a PointDataTree.
Definition: PointCountImpl.h:18
SharedPtr< const TreeBase > ConstPtr
Definition: Tree.h:41
const LeafNodeT * probeConstLeaf(const Coord &xyz) const
Return a pointer to the leaf node that contains the voxel coordinate xyz. If no LeafNode exists...
Definition: ValueAccessor.h:838
Definition: Exceptions.h:13
SharedPtr< const GridBase > ConstPtr
Definition: Grid.h:81
ValueAccessors are designed to help accelerate accesses into the OpenVDB Tree structures by storing c...
GridType
List of types that are currently supported by NanoVDB.
Definition: NanoVDB.h:219
uint32_t Index32
Definition: Types.h:52
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: ValueAccessor.h:455
bool diagonalizeSymmetricMatrix(const Mat3< T > &input, Mat3< T > &Q, Vec3< T > &D, unsigned int MAX_ITERATIONS=250)
Use Jacobi iterations to decompose a symmetric 3x3 matrix (diagonalize and compute eigenvectors) ...
Definition: Mat3.h:737
Vec4< int32_t > Vec4i
Definition: Vec4.h:559
SharedPtr< TreeBase > Ptr
Definition: Tree.h:40
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:121
Vec3< int32_t > Vec3i
Definition: Vec3.h:662
T dot(const Vec3< T > &v) const
Dot product.
Definition: Vec3.h:192
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:28
void setValue(const Coord &xyz, const ValueType &value)
Set a particular value at the given coordinate and mark the coordinate as active. ...
Definition: ValueAccessor.h:550
bool isValueOn(const Coord &xyz) const
Return the active state of the voxel at the given coordinates.
Definition: ValueAccessor.h:480
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:218
Gradient operators defined in index space of various orders.
Definition: Operators.h:99
Abstract base class for typed grids.
Definition: Grid.h:77