10 #ifndef OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
11 #define OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
19 #include <tbb/blocked_range.h>
20 #include <tbb/parallel_for.h>
21 #include <tbb/parallel_reduce.h>
22 #include <tbb/task_arena.h>
29 #include <type_traits>
53 template<
typename Gr
idType>
57 std::vector<Vec3s>&
points,
58 std::vector<Vec4I>&
quads,
59 double isovalue = 0.0);
82 template<
typename Gr
idType>
86 std::vector<Vec3s>&
points,
87 std::vector<Vec3I>& triangles,
88 std::vector<Vec4I>&
quads,
89 double isovalue = 0.0,
90 double adaptivity = 0.0,
91 bool relaxDisorientedTriangles =
true);
120 const size_t&
numQuads()
const {
return mNumQuads; }
144 inline bool trimQuads(
const size_t n,
bool reallocate =
false);
145 inline bool trimTrinagles(
const size_t n,
bool reallocate =
false);
151 size_t mNumQuads, mNumTriangles;
152 std::unique_ptr<openvdb::Vec4I[]> mQuads;
153 std::unique_ptr<openvdb::Vec3I[]> mTriangles;
154 std::unique_ptr<char[]> mQuadFlags, mTriangleFlags;
176 VolumeToMesh(
double isovalue = 0,
double adaptivity = 0,
bool relaxDisorientedTriangles =
true);
192 const std::vector<uint8_t>&
pointFlags()
const {
return mPointFlags; }
201 template<
typename InputGr
idType>
255 size_t mPointListSize, mSeamPointListSize, mPolygonPoolListSize;
256 double mIsovalue, mPrimAdaptivity, mSecAdaptivity;
263 bool mInvertSurfaceMask, mRelaxDisorientedTriangles;
265 std::unique_ptr<uint32_t[]> mQuantizedSeamPoints;
266 std::vector<uint8_t> mPointFlags;
280 const std::vector<Vec3d>&
points,
281 const std::vector<Vec3d>& normals)
287 if (points.empty())
return avgPos;
289 for (
size_t n = 0,
N = points.size();
n <
N; ++
n) {
293 avgPos /= double(points.size());
297 double m00=0,m01=0,m02=0,
304 for (
size_t n = 0, N = points.size();
n <
N; ++
n) {
306 const Vec3d& n_ref = normals[
n];
309 m00 += n_ref[0] * n_ref[0];
310 m11 += n_ref[1] * n_ref[1];
311 m22 += n_ref[2] * n_ref[2];
313 m01 += n_ref[0] * n_ref[1];
314 m02 += n_ref[0] * n_ref[2];
315 m12 += n_ref[1] * n_ref[2];
318 rhs += n_ref * n_ref.
dot(points[
n] - avgPos);
343 Mat3d D = Mat3d::identity();
351 for (
int i = 0; i < 3; ++i ) {
352 if (
std::abs(eigenValues[i]) < tolerance) {
356 D[i][i] = 1.0 / eigenValues[i];
363 return avgPos + pseudoInv * rhs;
378 namespace volume_to_mesh_internal {
380 template<
typename ValueType>
385 void operator()(
const tbb::blocked_range<size_t>&
range)
const {
387 for (
size_t n = range.begin(),
N = range.
end();
n <
N; ++
n) {
397 template<
typename ValueType>
401 const auto grainSize = std::max<size_t>(
402 length / tbb::this_task_arena::max_concurrency(), 1024);
403 const tbb::blocked_range<size_t>
range(0, length, grainSize);
404 tbb::parallel_for(range, FillArray<ValueType>(array, val), tbb::simple_partitioner());
409 enum { SIGNS = 0xFF,
EDGES = 0xE00,
INSIDE = 0x100,
410 XEDGE = 0x200, YEDGE = 0x400, ZEDGE = 0x800,
SEAM = 0x1000};
414 const bool sAdaptable[256] = {
415 1,1,1,1,1,0,1,1,1,1,0,1,1,1,1,1,1,1,0,1,0,0,0,1,0,1,0,1,0,1,0,1,
416 1,0,1,1,0,0,1,1,0,0,0,1,0,0,1,1,1,1,1,1,0,0,1,1,0,1,0,1,0,0,0,1,
417 1,0,0,0,1,0,1,1,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
418 1,0,1,1,1,0,1,1,0,0,0,0,1,0,1,1,1,1,1,1,1,0,1,1,0,0,0,0,0,0,0,1,
419 1,0,0,0,0,0,0,0,1,1,0,1,1,1,1,1,1,1,0,1,0,0,0,0,1,1,0,1,1,1,0,1,
420 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,1,1,0,1,0,0,0,1,
421 1,0,0,0,1,0,1,0,1,1,0,0,1,1,1,1,1,1,0,0,1,0,0,0,1,1,0,0,1,1,0,1,
422 1,0,1,0,1,0,1,0,1,0,0,0,1,0,1,1,1,1,1,1,1,0,1,1,1,1,0,1,1,1,1,1};
426 const unsigned char sAmbiguousFace[256] = {
427 0,0,0,0,0,5,0,0,0,0,5,0,0,0,0,0,0,0,1,0,0,5,1,0,4,0,0,0,4,0,0,0,
428 0,1,0,0,2,0,0,0,0,1,5,0,2,0,0,0,0,0,0,0,2,0,0,0,4,0,0,0,0,0,0,0,
429 0,0,2,2,0,5,0,0,3,3,0,0,0,0,0,0,6,6,0,0,6,0,0,0,0,0,0,0,0,0,0,0,
430 0,1,0,0,0,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
431 0,4,0,4,3,0,3,0,0,0,5,0,0,0,0,0,0,0,1,0,3,0,0,0,0,0,0,0,0,0,0,0,
432 6,0,6,0,0,0,0,0,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
433 0,4,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
434 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
440 const unsigned char sEdgeGroupTable[256][13] = {
441 {0,0,0,0,0,0,0,0,0,0,0,0,0},{1,1,0,0,1,0,0,0,0,1,0,0,0},{1,1,1,0,0,0,0,0,0,0,1,0,0},
442 {1,0,1,0,1,0,0,0,0,1,1,0,0},{1,0,1,1,0,0,0,0,0,0,0,1,0},{1,1,1,1,1,0,0,0,0,1,0,1,0},
443 {1,1,0,1,0,0,0,0,0,0,1,1,0},{1,0,0,1,1,0,0,0,0,1,1,1,0},{1,0,0,1,1,0,0,0,0,0,0,0,1},
444 {1,1,0,1,0,0,0,0,0,1,0,0,1},{1,1,1,1,1,0,0,0,0,0,1,0,1},{1,0,1,1,0,0,0,0,0,1,1,0,1},
445 {1,0,1,0,1,0,0,0,0,0,0,1,1},{1,1,1,0,0,0,0,0,0,1,0,1,1},{1,1,0,0,1,0,0,0,0,0,1,1,1},
446 {1,0,0,0,0,0,0,0,0,1,1,1,1},{1,0,0,0,0,1,0,0,1,1,0,0,0},{1,1,0,0,1,1,0,0,1,0,0,0,0},
447 {1,1,1,0,0,1,0,0,1,1,1,0,0},{1,0,1,0,1,1,0,0,1,0,1,0,0},{2,0,1,1,0,2,0,0,2,2,0,1,0},
448 {1,1,1,1,1,1,0,0,1,0,0,1,0},{1,1,0,1,0,1,0,0,1,1,1,1,0},{1,0,0,1,1,1,0,0,1,0,1,1,0},
449 {1,0,0,1,1,1,0,0,1,1,0,0,1},{1,1,0,1,0,1,0,0,1,0,0,0,1},{2,2,1,1,2,1,0,0,1,2,1,0,1},
450 {1,0,1,1,0,1,0,0,1,0,1,0,1},{1,0,1,0,1,1,0,0,1,1,0,1,1},{1,1,1,0,0,1,0,0,1,0,0,1,1},
451 {2,1,0,0,1,2,0,0,2,1,2,2,2},{1,0,0,0,0,1,0,0,1,0,1,1,1},{1,0,0,0,0,1,1,0,0,0,1,0,0},
452 {1,1,0,0,1,1,1,0,0,1,1,0,0},{1,1,1,0,0,1,1,0,0,0,0,0,0},{1,0,1,0,1,1,1,0,0,1,0,0,0},
453 {1,0,1,1,0,1,1,0,0,0,1,1,0},{2,2,2,1,1,1,1,0,0,1,2,1,0},{1,1,0,1,0,1,1,0,0,0,0,1,0},
454 {1,0,0,1,1,1,1,0,0,1,0,1,0},{2,0,0,2,2,1,1,0,0,0,1,0,2},{1,1,0,1,0,1,1,0,0,1,1,0,1},
455 {1,1,1,1,1,1,1,0,0,0,0,0,1},{1,0,1,1,0,1,1,0,0,1,0,0,1},{1,0,1,0,1,1,1,0,0,0,1,1,1},
456 {2,1,1,0,0,2,2,0,0,2,1,2,2},{1,1,0,0,1,1,1,0,0,0,0,1,1},{1,0,0,0,0,1,1,0,0,1,0,1,1},
457 {1,0,0,0,0,0,1,0,1,1,1,0,0},{1,1,0,0,1,0,1,0,1,0,1,0,0},{1,1,1,0,0,0,1,0,1,1,0,0,0},
458 {1,0,1,0,1,0,1,0,1,0,0,0,0},{1,0,1,1,0,0,1,0,1,1,1,1,0},{2,1,1,2,2,0,2,0,2,0,1,2,0},
459 {1,1,0,1,0,0,1,0,1,1,0,1,0},{1,0,0,1,1,0,1,0,1,0,0,1,0},{1,0,0,1,1,0,1,0,1,1,1,0,1},
460 {1,1,0,1,0,0,1,0,1,0,1,0,1},{2,1,2,2,1,0,2,0,2,1,0,0,2},{1,0,1,1,0,0,1,0,1,0,0,0,1},
461 {2,0,2,0,2,0,1,0,1,2,2,1,1},{2,2,2,0,0,0,1,0,1,0,2,1,1},{2,2,0,0,2,0,1,0,1,2,0,1,1},
462 {1,0,0,0,0,0,1,0,1,0,0,1,1},{1,0,0,0,0,0,1,1,0,0,0,1,0},{2,1,0,0,1,0,2,2,0,1,0,2,0},
463 {1,1,1,0,0,0,1,1,0,0,1,1,0},{1,0,1,0,1,0,1,1,0,1,1,1,0},{1,0,1,1,0,0,1,1,0,0,0,0,0},
464 {1,1,1,1,1,0,1,1,0,1,0,0,0},{1,1,0,1,0,0,1,1,0,0,1,0,0},{1,0,0,1,1,0,1,1,0,1,1,0,0},
465 {1,0,0,1,1,0,1,1,0,0,0,1,1},{1,1,0,1,0,0,1,1,0,1,0,1,1},{2,1,2,2,1,0,1,1,0,0,1,2,1},
466 {2,0,1,1,0,0,2,2,0,2,2,1,2},{1,0,1,0,1,0,1,1,0,0,0,0,1},{1,1,1,0,0,0,1,1,0,1,0,0,1},
467 {1,1,0,0,1,0,1,1,0,0,1,0,1},{1,0,0,0,0,0,1,1,0,1,1,0,1},{1,0,0,0,0,1,1,1,1,1,0,1,0},
468 {1,1,0,0,1,1,1,1,1,0,0,1,0},{2,1,1,0,0,2,2,1,1,1,2,1,0},{2,0,2,0,2,1,1,2,2,0,1,2,0},
469 {1,0,1,1,0,1,1,1,1,1,0,0,0},{2,2,2,1,1,2,2,1,1,0,0,0,0},{2,2,0,2,0,1,1,2,2,2,1,0,0},
470 {2,0,0,1,1,2,2,1,1,0,2,0,0},{2,0,0,1,1,1,1,2,2,1,0,1,2},{2,2,0,2,0,2,2,1,1,0,0,2,1},
471 {4,3,2,2,3,4,4,1,1,3,4,2,1},{3,0,2,2,0,1,1,3,3,0,1,2,3},{2,0,2,0,2,2,2,1,1,2,0,0,1},
472 {2,1,1,0,0,1,1,2,2,0,0,0,2},{3,1,0,0,1,2,2,3,3,1,2,0,3},{2,0,0,0,0,1,1,2,2,0,1,0,2},
473 {1,0,0,0,0,1,0,1,0,0,1,1,0},{1,1,0,0,1,1,0,1,0,1,1,1,0},{1,1,1,0,0,1,0,1,0,0,0,1,0},
474 {1,0,1,0,1,1,0,1,0,1,0,1,0},{1,0,1,1,0,1,0,1,0,0,1,0,0},{2,1,1,2,2,2,0,2,0,2,1,0,0},
475 {1,1,0,1,0,1,0,1,0,0,0,0,0},{1,0,0,1,1,1,0,1,0,1,0,0,0},{1,0,0,1,1,1,0,1,0,0,1,1,1},
476 {2,2,0,2,0,1,0,1,0,1,2,2,1},{2,2,1,1,2,2,0,2,0,0,0,1,2},{2,0,2,2,0,1,0,1,0,1,0,2,1},
477 {1,0,1,0,1,1,0,1,0,0,1,0,1},{2,2,2,0,0,1,0,1,0,1,2,0,1},{1,1,0,0,1,1,0,1,0,0,0,0,1},
478 {1,0,0,0,0,1,0,1,0,1,0,0,1},{1,0,0,0,0,0,0,1,1,1,1,1,0},{1,1,0,0,1,0,0,1,1,0,1,1,0},
479 {1,1,1,0,0,0,0,1,1,1,0,1,0},{1,0,1,0,1,0,0,1,1,0,0,1,0},{1,0,1,1,0,0,0,1,1,1,1,0,0},
480 {2,2,2,1,1,0,0,1,1,0,2,0,0},{1,1,0,1,0,0,0,1,1,1,0,0,0},{1,0,0,1,1,0,0,1,1,0,0,0,0},
481 {2,0,0,2,2,0,0,1,1,2,2,2,1},{2,1,0,1,0,0,0,2,2,0,1,1,2},{3,2,1,1,2,0,0,3,3,2,0,1,3},
482 {2,0,1,1,0,0,0,2,2,0,0,1,2},{2,0,1,0,1,0,0,2,2,1,1,0,2},{2,1,1,0,0,0,0,2,2,0,1,0,2},
483 {2,1,0,0,1,0,0,2,2,1,0,0,2},{1,0,0,0,0,0,0,1,1,0,0,0,1},{1,0,0,0,0,0,0,1,1,0,0,0,1},
484 {1,1,0,0,1,0,0,1,1,1,0,0,1},{2,1,1,0,0,0,0,2,2,0,1,0,2},{1,0,1,0,1,0,0,1,1,1,1,0,1},
485 {1,0,1,1,0,0,0,1,1,0,0,1,1},{2,1,1,2,2,0,0,1,1,1,0,1,2},{1,1,0,1,0,0,0,1,1,0,1,1,1},
486 {2,0,0,1,1,0,0,2,2,2,2,2,1},{1,0,0,1,1,0,0,1,1,0,0,0,0},{1,1,0,1,0,0,0,1,1,1,0,0,0},
487 {1,1,1,1,1,0,0,1,1,0,1,0,0},{1,0,1,1,0,0,0,1,1,1,1,0,0},{1,0,1,0,1,0,0,1,1,0,0,1,0},
488 {1,1,1,0,0,0,0,1,1,1,0,1,0},{1,1,0,0,1,0,0,1,1,0,1,1,0},{1,0,0,0,0,0,0,1,1,1,1,1,0},
489 {1,0,0,0,0,1,0,1,0,1,0,0,1},{1,1,0,0,1,1,0,1,0,0,0,0,1},{1,1,1,0,0,1,0,1,0,1,1,0,1},
490 {1,0,1,0,1,1,0,1,0,0,1,0,1},{1,0,1,1,0,1,0,1,0,1,0,1,1},{2,2,2,1,1,2,0,2,0,0,0,2,1},
491 {2,1,0,1,0,2,0,2,0,1,2,2,1},{2,0,0,2,2,1,0,1,0,0,1,1,2},{1,0,0,1,1,1,0,1,0,1,0,0,0},
492 {1,1,0,1,0,1,0,1,0,0,0,0,0},{2,1,2,2,1,2,0,2,0,1,2,0,0},{1,0,1,1,0,1,0,1,0,0,1,0,0},
493 {1,0,1,0,1,1,0,1,0,1,0,1,0},{1,1,1,0,0,1,0,1,0,0,0,1,0},{2,2,0,0,2,1,0,1,0,2,1,1,0},
494 {1,0,0,0,0,1,0,1,0,0,1,1,0},{1,0,0,0,0,1,1,1,1,0,1,0,1},{2,1,0,0,1,2,1,1,2,2,1,0,1},
495 {1,1,1,0,0,1,1,1,1,0,0,0,1},{2,0,2,0,2,1,2,2,1,1,0,0,2},{2,0,1,1,0,1,2,2,1,0,1,2,1},
496 {4,1,1,3,3,2,4,4,2,2,1,4,3},{2,2,0,2,0,2,1,1,2,0,0,1,2},{3,0,0,1,1,2,3,3,2,2,0,3,1},
497 {1,0,0,1,1,1,1,1,1,0,1,0,0},{2,2,0,2,0,1,2,2,1,1,2,0,0},{2,2,1,1,2,2,1,1,2,0,0,0,0},
498 {2,0,1,1,0,2,1,1,2,2,0,0,0},{2,0,2,0,2,2,1,1,2,0,2,1,0},{3,1,1,0,0,3,2,2,3,3,1,2,0},
499 {2,1,0,0,1,1,2,2,1,0,0,2,0},{2,0,0,0,0,2,1,1,2,2,0,1,0},{1,0,0,0,0,0,1,1,0,1,1,0,1},
500 {1,1,0,0,1,0,1,1,0,0,1,0,1},{1,1,1,0,0,0,1,1,0,1,0,0,1},{1,0,1,0,1,0,1,1,0,0,0,0,1},
501 {2,0,2,2,0,0,1,1,0,2,2,1,2},{3,1,1,2,2,0,3,3,0,0,1,3,2},{2,1,0,1,0,0,2,2,0,1,0,2,1},
502 {2,0,0,1,1,0,2,2,0,0,0,2,1},{1,0,0,1,1,0,1,1,0,1,1,0,0},{1,1,0,1,0,0,1,1,0,0,1,0,0},
503 {2,2,1,1,2,0,1,1,0,2,0,0,0},{1,0,1,1,0,0,1,1,0,0,0,0,0},{2,0,1,0,1,0,2,2,0,1,1,2,0},
504 {2,1,1,0,0,0,2,2,0,0,1,2,0},{2,1,0,0,1,0,2,2,0,1,0,2,0},{1,0,0,0,0,0,1,1,0,0,0,1,0},
505 {1,0,0,0,0,0,1,0,1,0,0,1,1},{1,1,0,0,1,0,1,0,1,1,0,1,1},{1,1,1,0,0,0,1,0,1,0,1,1,1},
506 {2,0,2,0,2,0,1,0,1,1,1,2,2},{1,0,1,1,0,0,1,0,1,0,0,0,1},{2,2,2,1,1,0,2,0,2,2,0,0,1},
507 {1,1,0,1,0,0,1,0,1,0,1,0,1},{2,0,0,2,2,0,1,0,1,1,1,0,2},{1,0,0,1,1,0,1,0,1,0,0,1,0},
508 {1,1,0,1,0,0,1,0,1,1,0,1,0},{2,2,1,1,2,0,2,0,2,0,2,1,0},{2,0,2,2,0,0,1,0,1,1,1,2,0},
509 {1,0,1,0,1,0,1,0,1,0,0,0,0},{1,1,1,0,0,0,1,0,1,1,0,0,0},{1,1,0,0,1,0,1,0,1,0,1,0,0},
510 {1,0,0,0,0,0,1,0,1,1,1,0,0},{1,0,0,0,0,1,1,0,0,1,0,1,1},{1,1,0,0,1,1,1,0,0,0,0,1,1},
511 {2,2,2,0,0,1,1,0,0,2,1,2,2},{2,0,1,0,1,2,2,0,0,0,2,1,1},{1,0,1,1,0,1,1,0,0,1,0,0,1},
512 {2,1,1,2,2,1,1,0,0,0,0,0,2},{2,1,0,1,0,2,2,0,0,1,2,0,1},{2,0,0,2,2,1,1,0,0,0,1,0,2},
513 {1,0,0,1,1,1,1,0,0,1,0,1,0},{1,1,0,1,0,1,1,0,0,0,0,1,0},{3,1,2,2,1,3,3,0,0,1,3,2,0},
514 {2,0,1,1,0,2,2,0,0,0,2,1,0},{1,0,1,0,1,1,1,0,0,1,0,0,0},{1,1,1,0,0,1,1,0,0,0,0,0,0},
515 {2,2,0,0,2,1,1,0,0,2,1,0,0},{1,0,0,0,0,1,1,0,0,0,1,0,0},{1,0,0,0,0,1,0,0,1,0,1,1,1},
516 {2,2,0,0,2,1,0,0,1,1,2,2,2},{1,1,1,0,0,1,0,0,1,0,0,1,1},{2,0,1,0,1,2,0,0,2,2,0,1,1},
517 {1,0,1,1,0,1,0,0,1,0,1,0,1},{3,1,1,3,3,2,0,0,2,2,1,0,3},{1,1,0,1,0,1,0,0,1,0,0,0,1},
518 {2,0,0,2,2,1,0,0,1,1,0,0,2},{1,0,0,1,1,1,0,0,1,0,1,1,0},{2,1,0,1,0,2,0,0,2,2,1,1,0},
519 {2,1,2,2,1,1,0,0,1,0,0,2,0},{2,0,1,1,0,2,0,0,2,2,0,1,0},{1,0,1,0,1,1,0,0,1,0,1,0,0},
520 {2,1,1,0,0,2,0,0,2,2,1,0,0},{1,1,0,0,1,1,0,0,1,0,0,0,0},{1,0,0,0,0,1,0,0,1,1,0,0,0},
521 {1,0,0,0,0,0,0,0,0,1,1,1,1},{1,1,0,0,1,0,0,0,0,0,1,1,1},{1,1,1,0,0,0,0,0,0,1,0,1,1},
522 {1,0,1,0,1,0,0,0,0,0,0,1,1},{1,0,1,1,0,0,0,0,0,1,1,0,1},{2,1,1,2,2,0,0,0,0,0,1,0,2},
523 {1,1,0,1,0,0,0,0,0,1,0,0,1},{1,0,0,1,1,0,0,0,0,0,0,0,1},{1,0,0,1,1,0,0,0,0,1,1,1,0},
524 {1,1,0,1,0,0,0,0,0,0,1,1,0},{2,1,2,2,1,0,0,0,0,1,0,2,0},{1,0,1,1,0,0,0,0,0,0,0,1,0},
525 {1,0,1,0,1,0,0,0,0,1,1,0,0},{1,1,1,0,0,0,0,0,0,0,1,0,0},{1,1,0,0,1,0,0,0,0,1,0,0,0},
526 {0,0,0,0,0,0,0,0,0,0,0,0,0}};
532 const Vec3d& p0,
const Vec3d& p1,
533 const Vec3d& p2,
const Vec3d& p3,
534 const double epsilon = 0.001)
539 const Vec3d centroid = (p0 + p1 + p2 + p3);
540 const double d = centroid.
dot(normal) * 0.25;
544 double absDist =
std::abs(p0.dot(normal) - d);
545 if (absDist > epsilon)
return false;
547 absDist =
std::abs(p1.dot(normal) - d);
548 if (absDist > epsilon)
return false;
550 absDist =
std::abs(p2.dot(normal) - d);
551 if (absDist > epsilon)
return false;
553 absDist =
std::abs(p3.dot(normal) - d);
554 if (absDist > epsilon)
return false;
567 MASK_FIRST_10_BITS = 0x000003FF,
568 MASK_DIRTY_BIT = 0x80000000,
569 MASK_INVALID_BIT = 0x40000000
573 packPoint(
const Vec3d& v)
578 assert(!(v.x() > 1.0) && !(v.y() > 1.0) && !(v.z() > 1.0));
579 assert(!(v.x() < 0.0) && !(v.y() < 0.0) && !(v.z() < 0.0));
581 data |= (uint32_t(v.x() * 1023.0) & MASK_FIRST_10_BITS) << 20;
582 data |= (uint32_t(v.y() * 1023.0) & MASK_FIRST_10_BITS) << 10;
583 data |= (uint32_t(v.z() * 1023.0) & MASK_FIRST_10_BITS);
589 unpackPoint(uint32_t data)
592 v.z() = double(data & MASK_FIRST_10_BITS) * 0.0009775171;
594 v.y() = double(data & MASK_FIRST_10_BITS) * 0.0009775171;
596 v.x() = double(data & MASK_FIRST_10_BITS) * 0.0009775171;
606 inline bool isBoolValue() {
return false; }
609 inline bool isBoolValue<bool>() {
return true; }
612 inline bool isInsideValue(T
value, T isovalue) {
return value < isovalue; }
615 inline bool isInsideValue<bool>(
bool value,
bool ) {
return value; }
620 template <
typename LeafT,
622 struct LeafBufferAccessor
625 LeafBufferAccessor(
const LeafT& leaf) : mData(leaf.
buffer().data()) {}
626 inline T get(
const Index idx)
const {
return mData[idx]; }
627 const T*
const mData;
630 template <
typename LeafT>
631 struct LeafBufferAccessor<LeafT, true>
634 LeafBufferAccessor(
const LeafT& leaf) : mLeaf(leaf) {}
635 inline T get(
const Index idx)
const {
return mLeaf.getValue(idx); }
641 template <
typename LeafT>
642 bool isInternalLeafCoord(
const Coord& ijk)
645 ijk[0] <
int(LeafT::DIM - 1) &&
646 ijk[1] <
int(LeafT::DIM - 1) &&
647 ijk[2] <
int(LeafT::DIM - 1);
652 template<
typename AccessorT,
typename ValueT>
654 getCellVertexValues(
const AccessorT& accessor,
656 std::array<ValueT, 8>&
values)
658 values[0] = ValueT(accessor.getValue(ijk));
660 values[1] = ValueT(accessor.getValue(ijk));
662 values[2] = ValueT(accessor.getValue(ijk));
664 values[3] = ValueT(accessor.getValue(ijk));
666 values[4] = ValueT(accessor.getValue(ijk));
668 values[5] = ValueT(accessor.getValue(ijk));
670 values[6] = ValueT(accessor.getValue(ijk));
672 values[7] = ValueT(accessor.getValue(ijk));
677 template<
typename LeafT,
typename ValueT>
679 getCellVertexValues(
const LeafT& leaf,
681 std::array<ValueT, 8>& values)
683 const LeafBufferAccessor<LeafT> acc(leaf);
685 values[0] = ValueT(acc.get(offset));
686 values[3] = ValueT(acc.get(offset + 1));
687 values[4] = ValueT(acc.get(offset + LeafT::DIM));
688 values[7] = ValueT(acc.get(offset + LeafT::DIM + 1));
689 values[1] = ValueT(acc.get(offset + (LeafT::DIM * LeafT::DIM)));
690 values[2] = ValueT(acc.get(offset + (LeafT::DIM * LeafT::DIM) + 1));
691 values[5] = ValueT(acc.get(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM));
692 values[6] = ValueT(acc.get(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM + 1));
696 template<
typename ValueType>
698 computeSignFlags(
const std::array<ValueType, 8>& values,
const ValueType iso)
701 signs |= isInsideValue(values[0], iso) ? 1u : 0u;
702 signs |= isInsideValue(values[1], iso) ? 2u : 0u;
703 signs |= isInsideValue(values[2], iso) ? 4u : 0u;
704 signs |= isInsideValue(values[3], iso) ? 8u : 0u;
705 signs |= isInsideValue(values[4], iso) ? 16u : 0u;
706 signs |= isInsideValue(values[5], iso) ? 32u : 0u;
707 signs |= isInsideValue(values[6], iso) ? 64u : 0u;
708 signs |= isInsideValue(values[7], iso) ? 128u : 0u;
709 return uint8_t(signs);
715 template<
typename AccessorT>
721 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 1u;
723 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 2u;
725 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 4u;
727 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 8u;
728 coord[1] += 1; coord[2] = ijk[2];
729 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 16u;
731 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 32u;
733 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 64u;
735 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 128u;
736 return uint8_t(signs);
742 template<
typename LeafT>
744 evalCellSigns(
const LeafT& leaf,
const Index offset,
typename LeafT::ValueType iso)
746 const LeafBufferAccessor<LeafT> acc(leaf);
749 if (isInsideValue(acc.get(offset), iso)) signs |= 1u;
750 if (isInsideValue(acc.get(offset + (LeafT::DIM * LeafT::DIM)), iso)) signs |= 2u;
751 if (isInsideValue(acc.get(offset + (LeafT::DIM * LeafT::DIM) + 1), iso)) signs |= 4u;
752 if (isInsideValue(acc.get(offset + 1), iso)) signs |= 8u;
753 if (isInsideValue(acc.get(offset + LeafT::DIM), iso)) signs |= 16u;
754 if (isInsideValue(acc.get(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM), iso)) signs |= 32u;
755 if (isInsideValue(acc.get(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM + 1), iso)) signs |= 64u;
756 if (isInsideValue(acc.get(offset + LeafT::DIM + 1), iso)) signs |= 128u;
757 return uint8_t(signs);
763 template<
class AccessorT>
765 correctCellSigns(uint8_t& signs,
767 const AccessorT& acc,
774 if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 3) signs = uint8_t(~signs);
778 if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 4) signs = uint8_t(~signs);
782 if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 1) signs = uint8_t(~signs);
786 if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 2) signs = uint8_t(~signs);
790 if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 6) signs = uint8_t(~signs);
794 if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 5) signs = uint8_t(~signs);
802 template<
class AccessorT>
804 isNonManifold(
const AccessorT& accessor,
const Coord& ijk,
807 const int hDim = dim >> 1;
811 p[0] = isInsideValue(accessor.getValue(coord), isovalue);
813 p[1] = isInsideValue(accessor.getValue(coord), isovalue);
815 p[2] = isInsideValue(accessor.getValue(coord), isovalue);
817 p[3] = isInsideValue(accessor.getValue(coord), isovalue);
818 coord[1] += dim; coord[2] = ijk[2];
819 p[4] = isInsideValue(accessor.getValue(coord), isovalue);
821 p[5] = isInsideValue(accessor.getValue(coord), isovalue);
823 p[6] = isInsideValue(accessor.getValue(coord), isovalue);
825 p[7] = isInsideValue(accessor.getValue(coord), isovalue);
829 if (p[0]) signs |= 1u;
830 if (p[1]) signs |= 2u;
831 if (p[2]) signs |= 4u;
832 if (p[3]) signs |= 8u;
833 if (p[4]) signs |= 16u;
834 if (p[5]) signs |= 32u;
835 if (p[6]) signs |= 64u;
836 if (p[7]) signs |= 128u;
837 if (!sAdaptable[signs])
return true;
842 const int i = ijk[0], ip = ijk[0] + hDim, ipp = ijk[0] + dim;
843 const int j = ijk[1], jp = ijk[1] + hDim, jpp = ijk[1] + dim;
844 const int k = ijk[2], kp = ijk[2] + hDim, kpp = ijk[2] + dim;
847 coord.reset(ip, j, k);
848 m = isInsideValue(accessor.getValue(coord), isovalue);
849 if (p[0] != m && p[1] != m)
return true;
852 coord.reset(ipp, j, kp);
853 m = isInsideValue(accessor.getValue(coord), isovalue);
854 if (p[1] != m && p[2] != m)
return true;
857 coord.reset(ip, j, kpp);
858 m = isInsideValue(accessor.getValue(coord), isovalue);
859 if (p[2] != m && p[3] != m)
return true;
862 coord.reset(i, j, kp);
863 m = isInsideValue(accessor.getValue(coord), isovalue);
864 if (p[0] != m && p[3] != m)
return true;
867 coord.reset(ip, jpp, k);
868 m = isInsideValue(accessor.getValue(coord), isovalue);
869 if (p[4] != m && p[5] != m)
return true;
872 coord.reset(ipp, jpp, kp);
873 m = isInsideValue(accessor.getValue(coord), isovalue);
874 if (p[5] != m && p[6] != m)
return true;
877 coord.reset(ip, jpp, kpp);
878 m = isInsideValue(accessor.getValue(coord), isovalue);
879 if (p[6] != m && p[7] != m)
return true;
882 coord.reset(i, jpp, kp);
883 m = isInsideValue(accessor.getValue(coord), isovalue);
884 if (p[7] != m && p[4] != m)
return true;
887 coord.reset(i, jp, k);
888 m = isInsideValue(accessor.getValue(coord), isovalue);
889 if (p[0] != m && p[4] != m)
return true;
892 coord.reset(ipp, jp, k);
893 m = isInsideValue(accessor.getValue(coord), isovalue);
894 if (p[1] != m && p[5] != m)
return true;
897 coord.reset(ipp, jp, kpp);
898 m = isInsideValue(accessor.getValue(coord), isovalue);
899 if (p[2] != m && p[6] != m)
return true;
903 coord.reset(i, jp, kpp);
904 m = isInsideValue(accessor.getValue(coord), isovalue);
905 if (p[3] != m && p[7] != m)
return true;
911 coord.reset(ip, jp, k);
912 m = isInsideValue(accessor.getValue(coord), isovalue);
913 if (p[0] != m && p[1] != m && p[4] != m && p[5] != m)
return true;
916 coord.reset(ipp, jp, kp);
917 m = isInsideValue(accessor.getValue(coord), isovalue);
918 if (p[1] != m && p[2] != m && p[5] != m && p[6] != m)
return true;
921 coord.reset(ip, jp, kpp);
922 m = isInsideValue(accessor.getValue(coord), isovalue);
923 if (p[2] != m && p[3] != m && p[6] != m && p[7] != m)
return true;
926 coord.reset(i, jp, kp);
927 m = isInsideValue(accessor.getValue(coord), isovalue);
928 if (p[0] != m && p[3] != m && p[4] != m && p[7] != m)
return true;
931 coord.reset(ip, j, kp);
932 m = isInsideValue(accessor.getValue(coord), isovalue);
933 if (p[0] != m && p[1] != m && p[2] != m && p[3] != m)
return true;
936 coord.reset(ip, jpp, kp);
937 m = isInsideValue(accessor.getValue(coord), isovalue);
938 if (p[4] != m && p[5] != m && p[6] != m && p[7] != m)
return true;
941 coord.reset(ip, jp, kp);
942 m = isInsideValue(accessor.getValue(coord), isovalue);
943 if (p[0] != m && p[1] != m && p[2] != m && p[3] != m &&
944 p[4] != m && p[5] != m && p[6] != m && p[7] != m)
return true;
953 template <
class LeafType>
955 mergeVoxels(LeafType& leaf,
const Coord&
start,
const int dim,
const int regionId)
958 const Coord
end = start.offsetBy(dim);
960 for (ijk[0] = start[0]; ijk[0] < end[0]; ++ijk[0]) {
961 for (ijk[1] = start[1]; ijk[1] < end[1]; ++ijk[1]) {
962 for (ijk[2] = start[2]; ijk[2] < end[2]; ++ijk[2]) {
963 leaf.setValueOnly(ijk, regionId);
972 template <
class LeafType>
974 isMergeable(
const LeafType& leaf,
979 if (adaptivity < 1e-6)
return false;
983 const Coord end = start.offsetBy(dim);
985 std::vector<VecT> norms;
986 for (ijk[0] = start[0]; ijk[0] < end[0]; ++ijk[0]) {
987 for (ijk[1] = start[1]; ijk[1] < end[1]; ++ijk[1]) {
988 for (ijk[2] = start[2]; ijk[2] < end[2]; ++ijk[2]) {
989 if (!leaf.isValueOn(ijk))
continue;
990 norms.push_back(leaf.getValue(ijk));
995 const size_t N = norms.size();
996 for (
size_t ni = 0; ni <
N; ++ni) {
997 VecT n_i = norms[ni];
998 for (
size_t nj = 0; nj <
N; ++nj) {
999 VecT n_j = norms[nj];
1000 if ((1.0 - n_i.dot(n_j)) > adaptivity)
return false;
1011 inline double evalZeroCrossing(
double v0,
double v1,
double iso) {
return (iso - v0) / (v1 -
v0); }
1016 computePoint(
const std::array<double, 8>& values,
1017 const unsigned char signs,
1018 const unsigned char edgeGroup,
1021 Vec3d avg(0.0, 0.0, 0.0);
1024 if (sEdgeGroupTable[signs][1] == edgeGroup) {
1025 avg[0] += evalZeroCrossing(values[0], values[1], iso);
1029 if (sEdgeGroupTable[signs][2] == edgeGroup) {
1031 avg[2] += evalZeroCrossing(values[1], values[2], iso);
1035 if (sEdgeGroupTable[signs][3] == edgeGroup) {
1036 avg[0] += evalZeroCrossing(values[3], values[2], iso);
1041 if (sEdgeGroupTable[signs][4] == edgeGroup) {
1042 avg[2] += evalZeroCrossing(values[0], values[3], iso);
1046 if (sEdgeGroupTable[signs][5] == edgeGroup) {
1047 avg[0] += evalZeroCrossing(values[4], values[5], iso);
1052 if (sEdgeGroupTable[signs][6] == edgeGroup) {
1055 avg[2] += evalZeroCrossing(values[5], values[6], iso);
1059 if (sEdgeGroupTable[signs][7] == edgeGroup) {
1060 avg[0] += evalZeroCrossing(values[7], values[6], iso);
1066 if (sEdgeGroupTable[signs][8] == edgeGroup) {
1068 avg[2] += evalZeroCrossing(values[4], values[7], iso);
1072 if (sEdgeGroupTable[signs][9] == edgeGroup) {
1073 avg[1] += evalZeroCrossing(values[0], values[4], iso);
1077 if (sEdgeGroupTable[signs][10] == edgeGroup) {
1079 avg[1] += evalZeroCrossing(values[1], values[5], iso);
1083 if (sEdgeGroupTable[signs][11] == edgeGroup) {
1085 avg[1] += evalZeroCrossing(values[2], values[6], iso);
1090 if (sEdgeGroupTable[signs][12] == edgeGroup) {
1091 avg[1] += evalZeroCrossing(values[3], values[7], iso);
1097 const double w = 1.0 / double(samples);
1108 computeMaskedPoint(Vec3d& avg,
1109 const std::array<double, 8>& values,
1110 const unsigned char signs,
1111 const unsigned char signsMask,
1112 const unsigned char edgeGroup,
1115 avg =
Vec3d(0.0, 0.0, 0.0);
1118 if (sEdgeGroupTable[signs][1] == edgeGroup &&
1119 sEdgeGroupTable[signsMask][1] == 0) {
1120 avg[0] += evalZeroCrossing(values[0], values[1], iso);
1124 if (sEdgeGroupTable[signs][2] == edgeGroup &&
1125 sEdgeGroupTable[signsMask][2] == 0) {
1127 avg[2] += evalZeroCrossing(values[1], values[2], iso);
1131 if (sEdgeGroupTable[signs][3] == edgeGroup &&
1132 sEdgeGroupTable[signsMask][3] == 0) {
1133 avg[0] += evalZeroCrossing(values[3], values[2], iso);
1138 if (sEdgeGroupTable[signs][4] == edgeGroup &&
1139 sEdgeGroupTable[signsMask][4] == 0) {
1140 avg[2] += evalZeroCrossing(values[0], values[3], iso);
1144 if (sEdgeGroupTable[signs][5] == edgeGroup &&
1145 sEdgeGroupTable[signsMask][5] == 0) {
1146 avg[0] += evalZeroCrossing(values[4], values[5], iso);
1151 if (sEdgeGroupTable[signs][6] == edgeGroup &&
1152 sEdgeGroupTable[signsMask][6] == 0) {
1155 avg[2] += evalZeroCrossing(values[5], values[6], iso);
1159 if (sEdgeGroupTable[signs][7] == edgeGroup &&
1160 sEdgeGroupTable[signsMask][7] == 0) {
1161 avg[0] += evalZeroCrossing(values[7], values[6], iso);
1167 if (sEdgeGroupTable[signs][8] == edgeGroup &&
1168 sEdgeGroupTable[signsMask][8] == 0) {
1170 avg[2] += evalZeroCrossing(values[4], values[7], iso);
1174 if (sEdgeGroupTable[signs][9] == edgeGroup &&
1175 sEdgeGroupTable[signsMask][9] == 0) {
1176 avg[1] += evalZeroCrossing(values[0], values[4], iso);
1180 if (sEdgeGroupTable[signs][10] == edgeGroup &&
1181 sEdgeGroupTable[signsMask][10] == 0) {
1183 avg[1] += evalZeroCrossing(values[1], values[5], iso);
1187 if (sEdgeGroupTable[signs][11] == edgeGroup &&
1188 sEdgeGroupTable[signsMask][11] == 0) {
1190 avg[1] += evalZeroCrossing(values[2], values[6], iso);
1195 if (sEdgeGroupTable[signs][12] == edgeGroup &&
1196 sEdgeGroupTable[signsMask][12] == 0) {
1197 avg[1] += evalZeroCrossing(values[3], values[7], iso);
1203 const double w = 1.0 / double(samples);
1214 computeWeightedPoint(
const Vec3d& p,
1215 const std::array<double, 8>& values,
1216 const unsigned char signs,
1217 const unsigned char edgeGroup,
1223 Vec3d avg(0.0, 0.0, 0.0);
1225 if (sEdgeGroupTable[signs][1] == edgeGroup) {
1226 avg[0] = evalZeroCrossing(values[0], values[1], iso);
1230 samples.push_back(avg);
1233 if (sEdgeGroupTable[signs][2] == edgeGroup) {
1236 avg[2] = evalZeroCrossing(values[1], values[2], iso);
1238 samples.push_back(avg);
1241 if (sEdgeGroupTable[signs][3] == edgeGroup) {
1242 avg[0] = evalZeroCrossing(values[3], values[2], iso);
1246 samples.push_back(avg);
1249 if (sEdgeGroupTable[signs][4] == edgeGroup) {
1252 avg[2] = evalZeroCrossing(values[0], values[3], iso);
1254 samples.push_back(avg);
1257 if (sEdgeGroupTable[signs][5] == edgeGroup) {
1258 avg[0] = evalZeroCrossing(values[4], values[5], iso);
1262 samples.push_back(avg);
1265 if (sEdgeGroupTable[signs][6] == edgeGroup) {
1268 avg[2] = evalZeroCrossing(values[5], values[6], iso);
1270 samples.push_back(avg);
1273 if (sEdgeGroupTable[signs][7] == edgeGroup) {
1274 avg[0] = evalZeroCrossing(values[7], values[6], iso);
1278 samples.push_back(avg);
1281 if (sEdgeGroupTable[signs][8] == edgeGroup) {
1284 avg[2] = evalZeroCrossing(values[4], values[7], iso);
1286 samples.push_back(avg);
1289 if (sEdgeGroupTable[signs][9] == edgeGroup) {
1291 avg[1] = evalZeroCrossing(values[0], values[4], iso);
1294 samples.push_back(avg);
1297 if (sEdgeGroupTable[signs][10] == edgeGroup) {
1299 avg[1] = evalZeroCrossing(values[1], values[5], iso);
1302 samples.push_back(avg);
1305 if (sEdgeGroupTable[signs][11] == edgeGroup) {
1307 avg[1] = evalZeroCrossing(values[2], values[6], iso);
1310 samples.push_back(avg);
1313 if (sEdgeGroupTable[signs][12] == edgeGroup) {
1315 avg[1] = evalZeroCrossing(values[3], values[7], iso);
1318 samples.push_back(avg);
1321 assert(!samples.empty());
1322 if (samples.size() == 1) {
1323 return samples.front();
1326 std::vector<double> weights;
1327 weights.reserve(samples.size());
1329 for (
const Vec3d&
s : samples) {
1330 weights.emplace_back((
s-p).lengthSqr());
1333 double minWeight = weights.front();
1334 double maxWeight = weights.front();
1336 for (
size_t i = 1, I = weights.size(); i < I; ++i) {
1337 minWeight =
std::min(minWeight, weights[i]);
1338 maxWeight =
std::max(maxWeight, weights[i]);
1341 const double offset = maxWeight + minWeight * 0.1;
1342 for (
size_t i = 0, I = weights.size(); i < I; ++i) {
1343 weights[i] = offset - weights[i];
1346 double weightSum = 0.0;
1347 for (
size_t i = 0, I = weights.size(); i < I; ++i) {
1348 weightSum += weights[i];
1352 for (
size_t i = 0, I = samples.size(); i < I; ++i) {
1353 avg += samples[i] * (weights[i] / weightSum);
1363 computeCellPoints(std::array<Vec3d, 4>&
points,
1364 const std::array<double, 8>& values,
1365 const unsigned char signs,
1369 for (
size_t n = 1, N = sEdgeGroupTable[signs][0] + 1;
n <
N; ++
n, ++
offset) {
1371 points[
offset] = computePoint(values, signs, uint8_t(
n), iso);
1381 matchEdgeGroup(
unsigned char groupId,
unsigned char lhsSigns,
unsigned char rhsSigns)
1384 for (
size_t i = 1; i <= 12; ++i) {
1385 if (sEdgeGroupTable[lhsSigns][i] == groupId && sEdgeGroupTable[rhsSigns][i] != 0) {
1386 id = sEdgeGroupTable[rhsSigns][i];
1399 computeCellPoints(std::array<Vec3d, 4>& points,
1400 std::array<bool, 4>& weightedPointMask,
1401 const std::array<double, 8>& lhsValues,
1402 const std::array<double, 8>& rhsValues,
1403 const unsigned char lhsSigns,
1404 const unsigned char rhsSigns,
1406 const size_t pointIdx,
1407 const uint32_t * seamPointArray)
1410 for (
size_t n = 1, N = sEdgeGroupTable[lhsSigns][0] + 1;
n <
N; ++
n, ++
offset)
1413 const int id = matchEdgeGroup(uint8_t(
n), lhsSigns, rhsSigns);
1417 const unsigned char e = uint8_t(
id);
1418 const uint32_t quantizedPoint = seamPointArray[pointIdx + (
id - 1)];
1420 if ((quantizedPoint & MASK_DIRTY_BIT) && !(quantizedPoint & MASK_INVALID_BIT)) {
1421 const Vec3d p = unpackPoint(quantizedPoint);
1422 points[
offset] = computeWeightedPoint(p, rhsValues, rhsSigns, e, iso);
1423 weightedPointMask[
offset] =
true;
1425 points[
offset] = computePoint(rhsValues, rhsSigns, e, iso);
1426 weightedPointMask[
offset] =
false;
1430 points[
offset] = computePoint(lhsValues, lhsSigns, uint8_t(
n), iso);
1431 weightedPointMask[
offset] =
false;
1438 template <
typename InputTreeType>
1439 struct ComputePoints
1441 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
1445 using Int16LeafNodeType =
typename Int16TreeType::LeafNodeType;
1448 using Index32LeafNodeType =
typename Index32TreeType::LeafNodeType;
1450 ComputePoints(
Vec3s * pointArray,
1451 const InputTreeType& inputTree,
1452 const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
1453 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
1454 const std::unique_ptr<
Index32[]>& leafNodeOffsets,
1455 const math::Transform& xform,
1458 void setRefData(
const InputTreeType& refInputTree,
1459 const Index32TreeType& refPointIndexTree,
1460 const Int16TreeType& refSignFlagsTree,
1461 const uint32_t * quantizedSeamLinePoints,
1462 uint8_t * seamLinePointsFlags);
1464 void operator()(
const tbb::blocked_range<size_t>&)
const;
1467 Vec3s *
const mPoints;
1468 InputTreeType
const *
const mInputTree;
1469 Index32LeafNodeType *
const *
const mPointIndexNodes;
1470 Int16LeafNodeType
const *
const *
const mSignFlagsNodes;
1471 Index32 const *
const mNodeOffsets;
1472 math::Transform
const mTransform;
1473 double const mIsovalue;
1475 InputTreeType
const * mRefInputTree;
1476 Index32TreeType
const * mRefPointIndexTree;
1477 Int16TreeType
const * mRefSignFlagsTree;
1478 uint32_t
const * mQuantizedSeamLinePoints;
1479 uint8_t * mSeamLinePointsFlags;
1483 template <
typename InputTreeType>
1484 ComputePoints<InputTreeType>::ComputePoints(
1486 const InputTreeType& inputTree,
1487 const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
1488 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
1489 const std::unique_ptr<
Index32[]>& leafNodeOffsets,
1490 const math::Transform& xform,
1492 : mPoints(pointArray)
1493 , mInputTree(&inputTree)
1494 , mPointIndexNodes(pointIndexLeafNodes.data())
1495 , mSignFlagsNodes(signFlagsLeafNodes.data())
1496 , mNodeOffsets(leafNodeOffsets.
get())
1499 , mRefInputTree(nullptr)
1500 , mRefPointIndexTree(nullptr)
1501 , mRefSignFlagsTree(nullptr)
1502 , mQuantizedSeamLinePoints(nullptr)
1503 , mSeamLinePointsFlags(nullptr)
1507 template <
typename InputTreeType>
1509 ComputePoints<InputTreeType>::setRefData(
1510 const InputTreeType& refInputTree,
1511 const Index32TreeType& refPointIndexTree,
1512 const Int16TreeType& refSignFlagsTree,
1513 const uint32_t * quantizedSeamLinePoints,
1514 uint8_t * seamLinePointsFlags)
1516 mRefInputTree = &refInputTree;
1517 mRefPointIndexTree = &refPointIndexTree;
1518 mRefSignFlagsTree = &refSignFlagsTree;
1519 mQuantizedSeamLinePoints = quantizedSeamLinePoints;
1520 mSeamLinePointsFlags = seamLinePointsFlags;
1523 template <
typename InputTreeType>
1525 ComputePoints<InputTreeType>::operator()(
const tbb::blocked_range<size_t>& range)
const
1527 using InputTreeAccessor = tree::ValueAccessor<const InputTreeType>;
1528 using Index32TreeAccessor = tree::ValueAccessor<const Index32TreeType>;
1529 using Int16TreeAccessor = tree::ValueAccessor<const Int16TreeType>;
1534 using IndexArrayMap = std::map<IndexType, IndexArray>;
1536 InputTreeAccessor inputAcc(*mInputTree);
1540 std::array<Vec3d, 4>
points;
1541 std::array<bool, 4> weightedPointMask;
1542 std::array<double, 8>
values, refValues;
1543 const double iso = mIsovalue;
1547 std::unique_ptr<InputTreeAccessor> refInputAcc;
1548 std::unique_ptr<Index32TreeAccessor> refPointIndexAcc;
1549 std::unique_ptr<Int16TreeAccessor> refSignFlagsAcc;
1551 const bool hasReferenceData = mRefInputTree && mRefPointIndexTree && mRefSignFlagsTree;
1553 if (hasReferenceData) {
1554 refInputAcc.reset(
new InputTreeAccessor(*mRefInputTree));
1555 refPointIndexAcc.reset(
new Index32TreeAccessor(*mRefPointIndexTree));
1556 refSignFlagsAcc.reset(
new Int16TreeAccessor(*mRefSignFlagsTree));
1559 for (
size_t n = range.begin(), N = range.end();
n !=
N; ++
n)
1561 Index32LeafNodeType& pointIndexNode = *mPointIndexNodes[
n];
1562 const Coord& origin = pointIndexNode.origin();
1564 const Int16LeafNodeType& signFlagsNode = *mSignFlagsNodes[
n];
1565 const InputLeafNodeType * inputNode = inputAcc.probeConstLeaf(origin);
1568 const InputLeafNodeType * refInputNode =
nullptr;
1569 const Index32LeafNodeType * refPointIndexNode =
nullptr;
1570 const Int16LeafNodeType * refSignFlagsNode =
nullptr;
1572 if (hasReferenceData) {
1573 refInputNode = refInputAcc->probeConstLeaf(origin);
1574 refPointIndexNode = refPointIndexAcc->probeConstLeaf(origin);
1575 refSignFlagsNode = refSignFlagsAcc->probeConstLeaf(origin);
1579 IndexArrayMap regions;
1581 auto*
const pidxData = pointIndexNode.buffer().data();
1582 const auto*
const sfData = signFlagsNode.buffer().data();
1584 for (
auto it = pointIndexNode.beginValueOn(); it; ++it)
1586 const Index offset = it.pos();
1591 regions[
id].push_back(offset);
1599 const uint8_t signs = uint8_t(SIGNS & flags);
1600 uint8_t refSigns = 0;
1602 if ((flags & SEAM) && refPointIndexNode && refSignFlagsNode) {
1603 if (refSignFlagsNode->isValueOn(offset)) {
1604 refSigns = uint8_t(SIGNS & refSignFlagsNode->getValue(offset));
1608 ijk = Index32LeafNodeType::offsetToLocalCoord(offset);
1610 const bool inclusiveCell = inputNode && isInternalLeafCoord<InputLeafNodeType>(ijk);
1614 if (inclusiveCell) getCellVertexValues(*inputNode, offset, values);
1615 else getCellVertexValues(inputAcc, ijk, values);
1617 size_t count, weightcount;
1619 if (refSigns == 0) {
1620 count = computeCellPoints(points, values, signs, iso);
1623 if (inclusiveCell && refInputNode) {
1624 getCellVertexValues(*refInputNode, offset, refValues);
1626 getCellVertexValues(*refInputAcc, ijk, refValues);
1628 count = computeCellPoints(points, weightedPointMask, values, refValues, signs, refSigns,
1629 iso, refPointIndexNode->getValue(offset), mQuantizedSeamLinePoints);
1630 weightcount = count;
1633 xyz = ijk.asVec3d();
1635 for (
size_t i = 0; i < count; ++i) {
1637 Vec3d& point = points[i];
1640 if (!std::isfinite(point[0]) ||
1641 !std::isfinite(point[1]) ||
1642 !std::isfinite(point[2]))
1645 "VolumeToMesh encountered NaNs or infs in the input VDB!"
1646 " Hint: Check the input and consider using the \"Diagnostics\" tool "
1647 "to detect and resolve the NaNs.");
1651 point = mTransform.indexToWorld(point);
1653 Vec3s& pos = mPoints[pointOffset];
1654 pos[0] =
float(point[0]);
1655 pos[1] =
float(point[1]);
1656 pos[2] =
float(point[2]);
1658 if (mSeamLinePointsFlags && weightcount && weightedPointMask[i]) {
1659 mSeamLinePointsFlags[pointOffset] = uint8_t(1);
1667 for (
auto it = regions.begin(); it != regions.end(); ++it)
1672 for (
size_t i = 0, I = voxels.size(); i < I; ++i) {
1674 const Index offset = voxels[i];
1675 ijk = Index32LeafNodeType::offsetToLocalCoord(offset);
1677 const bool inclusiveCell = inputNode && isInternalLeafCoord<InputLeafNodeType>(ijk);
1681 pidxData[
offset] = pointOffset;
1683 const uint8_t signs = uint8_t(SIGNS & sfData[offset]);
1685 if (inclusiveCell) getCellVertexValues(*inputNode, offset, values);
1686 else getCellVertexValues(inputAcc, ijk, values);
1688 computeCellPoints(points, values, signs, iso);
1690 avg[0] += double(ijk[0]) + points[0][0];
1691 avg[1] += double(ijk[1]) + points[0][1];
1692 avg[2] += double(ijk[2]) + points[0][2];
1695 if (voxels.size() > 1) {
1696 const double w = 1.0 / double(voxels.size());
1700 avg = mTransform.indexToWorld(avg);
1702 Vec3s& pos = mPoints[pointOffset];
1703 pos[0] =
float(avg[0]);
1704 pos[1] =
float(avg[1]);
1705 pos[2] =
float(avg[2]);
1716 template <
typename InputTreeType>
1717 struct SeamLineWeights
1719 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
1723 using Int16LeafNodeType =
typename Int16TreeType::LeafNodeType;
1726 using Index32LeafNodeType =
typename Index32TreeType::LeafNodeType;
1728 SeamLineWeights(
const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
1729 const InputTreeType& inputTree,
1730 const Index32TreeType& refPointIndexTree,
1731 const Int16TreeType& refSignFlagsTree,
1732 uint32_t * quantizedPoints,
1734 : mSignFlagsNodes(signFlagsLeafNodes.data())
1735 , mInputTree(&inputTree)
1736 , mRefPointIndexTree(&refPointIndexTree)
1737 , mRefSignFlagsTree(&refSignFlagsTree)
1738 , mQuantizedPoints(quantizedPoints)
1743 void operator()(
const tbb::blocked_range<size_t>& range)
const
1745 tree::ValueAccessor<const InputTreeType> inputTreeAcc(*mInputTree);
1746 tree::ValueAccessor<const Index32TreeType> pointIndexTreeAcc(*mRefPointIndexTree);
1747 tree::ValueAccessor<const Int16TreeType> signFlagsTreeAcc(*mRefSignFlagsTree);
1749 std::array<double, 8>
values;
1750 const double iso = double(mIsovalue);
1754 for (
size_t n = range.begin(), N = range.end();
n !=
N; ++
n) {
1756 const Int16LeafNodeType& signFlagsNode = *mSignFlagsNodes[
n];
1757 const Coord& origin = signFlagsNode.origin();
1759 const Int16LeafNodeType * refSignNode = signFlagsTreeAcc.probeConstLeaf(origin);
1760 if (!refSignNode)
continue;
1762 const Index32LeafNodeType* refPointIndexNode =
1763 pointIndexTreeAcc.probeConstLeaf(origin);
1764 if (!refPointIndexNode)
continue;
1766 const InputLeafNodeType* inputNode = inputTreeAcc.probeConstLeaf(origin);
1768 const auto*
const sfData = signFlagsNode.buffer().data();
1769 const auto*
const rfIdxData = refPointIndexNode->buffer().data();
1770 const auto*
const rsData = refSignNode->buffer().data();
1772 for (
auto it = signFlagsNode.cbeginValueOn(); it; ++it)
1774 const Index offset = it.pos();
1777 ijk = Index32LeafNodeType::offsetToLocalCoord(offset);
1779 const bool inclusiveCell = inputNode && isInternalLeafCoord<InputLeafNodeType>(ijk);
1783 if ((flags & SEAM) && refSignNode->isValueOn(offset)) {
1785 const uint8_t lhsSigns = uint8_t(SIGNS & flags);
1786 const uint8_t rhsSigns = uint8_t(SIGNS & rsData[offset]);
1788 if (inclusiveCell) getCellVertexValues(*inputNode, offset, values);
1789 else getCellVertexValues(inputTreeAcc, ijk, values);
1791 for (
unsigned i = 1, I = sEdgeGroupTable[lhsSigns][0] + 1; i < I; ++i) {
1793 const int id = matchEdgeGroup(uint8_t(i), lhsSigns, rhsSigns);
1797 uint32_t& data = mQuantizedPoints[rfIdxData[
offset] + (
id - 1)];
1799 if (!(data & MASK_DIRTY_BIT)) {
1801 const int samples = computeMaskedPoint(
1802 pos, values, lhsSigns, rhsSigns, uint8_t(i), iso);
1804 if (samples > 0) data = packPoint(pos);
1805 else data = MASK_INVALID_BIT;
1807 data |= MASK_DIRTY_BIT;
1817 Int16LeafNodeType
const *
const *
const mSignFlagsNodes;
1818 InputTreeType
const *
const mInputTree;
1819 Index32TreeType
const *
const mRefPointIndexTree;
1820 Int16TreeType
const *
const mRefSignFlagsTree;
1821 uint32_t *
const mQuantizedPoints;
1822 InputValueType
const mIsovalue;
1826 template <
typename TreeType>
1827 struct SetSeamLineFlags
1829 using LeafNodeType =
typename TreeType::LeafNodeType;
1831 SetSeamLineFlags(
const std::vector<LeafNodeType*>& signFlagsLeafNodes,
1832 const TreeType& refSignFlagsTree)
1833 : mSignFlagsNodes(signFlagsLeafNodes.data())
1834 , mRefSignFlagsTree(&refSignFlagsTree)
1838 void operator()(
const tbb::blocked_range<size_t>& range)
const
1840 tree::ValueAccessor<const TreeType> refSignFlagsTreeAcc(*mRefSignFlagsTree);
1842 for (
size_t n = range.begin(), N = range.end();
n !=
N; ++
n) {
1844 LeafNodeType& signFlagsNode = *mSignFlagsNodes[
n];
1845 const Coord& origin = signFlagsNode.origin();
1847 const LeafNodeType* refSignNode = refSignFlagsTreeAcc.probeConstLeaf(origin);
1848 if (!refSignNode)
continue;
1850 const auto*
const rsData = refSignNode->buffer().data();
1851 auto*
const sfData = signFlagsNode.buffer().data();
1853 for (
auto it = signFlagsNode.cbeginValueOn(); it; ++it) {
1854 const Index offset = it.pos();
1856 const uint8_t rhsSigns = uint8_t(rsData[offset] & SIGNS);
1858 if (sEdgeGroupTable[rhsSigns][0] > 0) {
1861 const uint8_t lhsSigns = uint8_t(value & SIGNS);
1863 if (rhsSigns != lhsSigns) {
1874 LeafNodeType *
const *
const mSignFlagsNodes;
1875 TreeType
const *
const mRefSignFlagsTree;
1879 template <
typename BoolTreeType,
typename SignDataType>
1880 struct TransferSeamLineFlags
1882 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
1885 using SignDataLeafNodeType =
typename SignDataTreeType::LeafNodeType;
1887 TransferSeamLineFlags(
const std::vector<SignDataLeafNodeType*>& signFlagsLeafNodes,
1888 const BoolTreeType& maskTree)
1889 : mSignFlagsNodes(signFlagsLeafNodes.data())
1890 , mMaskTree(&maskTree)
1894 void operator()(
const tbb::blocked_range<size_t>& range)
const
1896 tree::ValueAccessor<const BoolTreeType> maskAcc(*mMaskTree);
1898 for (
size_t n = range.begin(), N = range.end();
n !=
N; ++
n) {
1900 SignDataLeafNodeType& signFlagsNode = *mSignFlagsNodes[
n];
1901 const Coord& origin = signFlagsNode.origin();
1903 const BoolLeafNodeType * maskNode = maskAcc.probeConstLeaf(origin);
1904 if (!maskNode)
continue;
1906 auto*
const sfData = signFlagsNode.buffer().data();
1908 for (
auto it = signFlagsNode.cbeginValueOn(); it; ++it) {
1909 const Index offset = it.pos();
1911 if (maskNode->isValueOn(offset)) {
1919 SignDataLeafNodeType *
const *
const mSignFlagsNodes;
1920 BoolTreeType
const *
const mMaskTree;
1924 template <
typename TreeType>
1925 struct MaskSeamLineVoxels
1927 using LeafNodeType =
typename TreeType::LeafNodeType;
1930 MaskSeamLineVoxels(
const std::vector<LeafNodeType*>& signFlagsLeafNodes,
1931 const TreeType& signFlagsTree,
1933 : mSignFlagsNodes(signFlagsLeafNodes.data())
1934 , mSignFlagsTree(&signFlagsTree)
1940 MaskSeamLineVoxels(MaskSeamLineVoxels& rhs,
tbb::split)
1941 : mSignFlagsNodes(rhs.mSignFlagsNodes)
1942 , mSignFlagsTree(rhs.mSignFlagsTree)
1948 void join(MaskSeamLineVoxels& rhs) {
mMask->merge(*rhs.mMask); }
1950 void operator()(
const tbb::blocked_range<size_t>& range)
1952 using ValueOnCIter =
typename LeafNodeType::ValueOnCIter;
1955 tree::ValueAccessor<const TreeType> signFlagsAcc(*mSignFlagsTree);
1956 tree::ValueAccessor<BoolTreeType> maskAcc(*
mMask);
1959 for (
size_t n = range.begin(), N = range.end();
n !=
N; ++
n) {
1961 LeafNodeType& signFlagsNode = *mSignFlagsNodes[
n];
1962 auto*
const sfData = signFlagsNode.buffer().data();
1964 for (ValueOnCIter it = signFlagsNode.cbeginValueOn(); it; ++it) {
1966 const ValueType flags = sfData[it.pos()];
1968 if (!(flags & SEAM) && (flags & EDGES)) {
1970 ijk = it.getCoord();
1972 bool isSeamLineVoxel =
false;
1974 if (flags & XEDGE) {
1976 isSeamLineVoxel = (signFlagsAcc.getValue(ijk) & SEAM);
1978 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
1980 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
1984 if (!isSeamLineVoxel && flags & YEDGE) {
1986 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
1988 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
1990 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
1994 if (!isSeamLineVoxel && flags & ZEDGE) {
1996 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
1998 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
2000 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
2004 if (isSeamLineVoxel) {
2005 maskAcc.setValue(it.getCoord(),
true);
2014 LeafNodeType *
const *
const mSignFlagsNodes;
2015 TreeType
const *
const mSignFlagsTree;
2016 BoolTreeType mTempMask;
2017 BoolTreeType *
const mMask;
2021 template<
typename SignDataTreeType>
2023 markSeamLineData(SignDataTreeType& signFlagsTree,
const SignDataTreeType& refSignFlagsTree)
2026 using SignDataLeafNodeType =
typename SignDataTreeType::LeafNodeType;
2029 std::vector<SignDataLeafNodeType*> signFlagsLeafNodes;
2030 signFlagsTree.getNodes(signFlagsLeafNodes);
2032 const tbb::blocked_range<size_t> nodeRange(0, signFlagsLeafNodes.size());
2035 SetSeamLineFlags<SignDataTreeType>(signFlagsLeafNodes, refSignFlagsTree));
2037 BoolTreeType seamLineMaskTree(
false);
2039 MaskSeamLineVoxels<SignDataTreeType>
2040 maskSeamLine(signFlagsLeafNodes, signFlagsTree, seamLineMaskTree);
2042 tbb::parallel_reduce(nodeRange, maskSeamLine);
2045 TransferSeamLineFlags<BoolTreeType, SignDataType>(signFlagsLeafNodes, seamLineMaskTree));
2052 template <
typename InputGr
idType>
2053 struct MergeVoxelRegions
2055 using InputTreeType =
typename InputGridType::TreeType;
2056 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
2060 using FloatLeafNodeType =
typename FloatTreeType::LeafNodeType;
2061 using FloatGridType = Grid<FloatTreeType>;
2064 using Int16LeafNodeType =
typename Int16TreeType::LeafNodeType;
2067 using Index32LeafNodeType =
typename Index32TreeType::LeafNodeType;
2069 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
2070 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
2073 using MaskLeafNodeType =
typename MaskTreeType::LeafNodeType;
2075 MergeVoxelRegions(
const InputGridType& inputGrid,
2076 const Index32TreeType& pointIndexTree,
2077 const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
2078 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
2081 bool invertSurfaceOrientation);
2083 void setSpatialAdaptivity(
const FloatGridType& grid)
2085 mSpatialAdaptivityTree = &grid.tree();
2086 mSpatialAdaptivityTransform = &grid.transform();
2089 void setAdaptivityMask(
const BoolTreeType&
mask)
2094 void setRefSignFlagsData(
const Int16TreeType& signFlagsData,
float internalAdaptivity)
2096 mRefSignFlagsTree = &signFlagsData;
2097 mInternalAdaptivity = internalAdaptivity;
2100 void operator()(
const tbb::blocked_range<size_t>&)
const;
2103 InputTreeType
const *
const mInputTree;
2104 math::Transform
const *
const mInputTransform;
2106 Index32TreeType
const *
const mPointIndexTree;
2107 Index32LeafNodeType *
const *
const mPointIndexNodes;
2108 Int16LeafNodeType
const *
const *
const mSignFlagsNodes;
2110 InputValueType mIsovalue;
2111 float mSurfaceAdaptivity, mInternalAdaptivity;
2112 bool mInvertSurfaceOrientation;
2114 FloatTreeType
const * mSpatialAdaptivityTree;
2115 BoolTreeType
const * mMaskTree;
2116 Int16TreeType
const * mRefSignFlagsTree;
2117 math::Transform
const * mSpatialAdaptivityTransform;
2121 template <
typename InputGr
idType>
2122 MergeVoxelRegions<InputGridType>::MergeVoxelRegions(
2123 const InputGridType& inputGrid,
2124 const Index32TreeType& pointIndexTree,
2125 const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
2126 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
2129 bool invertSurfaceOrientation)
2130 : mInputTree(&inputGrid.tree())
2131 , mInputTransform(&inputGrid.
transform())
2132 , mPointIndexTree(&pointIndexTree)
2133 , mPointIndexNodes(pointIndexLeafNodes.data())
2134 , mSignFlagsNodes(signFlagsLeafNodes.data())
2136 , mSurfaceAdaptivity(adaptivity)
2137 , mInternalAdaptivity(adaptivity)
2138 , mInvertSurfaceOrientation(invertSurfaceOrientation)
2139 , mSpatialAdaptivityTree(nullptr)
2140 , mMaskTree(nullptr)
2141 , mRefSignFlagsTree(nullptr)
2142 , mSpatialAdaptivityTransform(nullptr)
2147 template <
typename InputGr
idType>
2149 MergeVoxelRegions<InputGridType>::operator()(
const tbb::blocked_range<size_t>& range)
const
2151 using Vec3sType = math::Vec3<float>;
2154 using InputTreeAccessor = tree::ValueAccessor<const InputTreeType>;
2155 using FloatTreeAccessor = tree::ValueAccessor<const FloatTreeType>;
2156 using Index32TreeAccessor = tree::ValueAccessor<const Index32TreeType>;
2157 using Int16TreeAccessor = tree::ValueAccessor<const Int16TreeType>;
2158 using BoolTreeAccessor = tree::ValueAccessor<const BoolTreeType>;
2160 std::unique_ptr<FloatTreeAccessor> spatialAdaptivityAcc;
2161 if (mSpatialAdaptivityTree && mSpatialAdaptivityTransform) {
2162 spatialAdaptivityAcc.reset(
new FloatTreeAccessor(*mSpatialAdaptivityTree));
2165 std::unique_ptr<BoolTreeAccessor> maskAcc;
2167 maskAcc.reset(
new BoolTreeAccessor(*mMaskTree));
2170 std::unique_ptr<Int16TreeAccessor> refSignFlagsAcc;
2171 if (mRefSignFlagsTree) {
2172 refSignFlagsAcc.reset(
new Int16TreeAccessor(*mRefSignFlagsTree));
2175 InputTreeAccessor inputAcc(*mInputTree);
2176 Index32TreeAccessor pointIndexAcc(*mPointIndexTree);
2178 MaskLeafNodeType
mask;
2180 const bool invertGradientDir = mInvertSurfaceOrientation || isBoolValue<InputValueType>();
2181 std::unique_ptr<Vec3sLeafNodeType> gradientNode;
2184 const int LeafDim = InputLeafNodeType::DIM;
2186 for (
size_t n = range.begin(), N = range.end();
n !=
N; ++
n) {
2188 mask.setValuesOff();
2190 const Int16LeafNodeType& signFlagsNode = *mSignFlagsNodes[
n];
2191 Index32LeafNodeType& pointIndexNode = *mPointIndexNodes[
n];
2193 const Coord& origin = pointIndexNode.origin();
2194 end = origin.offsetBy(LeafDim);
2198 const BoolLeafNodeType* maskLeaf = maskAcc->probeConstLeaf(origin);
2199 if (maskLeaf !=
nullptr) {
2200 for (
auto it = maskLeaf->cbeginValueOn(); it; ++it)
2202 mask.setActiveState(it.getCoord() & ~1u,
true);
2207 float adaptivity = (refSignFlagsAcc && !refSignFlagsAcc->probeConstLeaf(origin)) ?
2208 mInternalAdaptivity : mSurfaceAdaptivity;
2210 bool useGradients = adaptivity < 1.0f;
2213 FloatLeafNodeType adaptivityLeaf(origin, adaptivity);
2215 if (spatialAdaptivityAcc) {
2216 useGradients =
false;
2217 for (
Index offset = 0; offset < FloatLeafNodeType::NUM_VALUES; ++
offset) {
2218 ijk = adaptivityLeaf.offsetToGlobalCoord(offset);
2219 ijk = mSpatialAdaptivityTransform->worldToIndexCellCentered(
2220 mInputTransform->indexToWorld(ijk));
2221 float weight = spatialAdaptivityAcc->getValue(ijk);
2222 float adaptivityValue = weight * adaptivity;
2223 if (adaptivityValue < 1.0
f) useGradients =
true;
2224 adaptivityLeaf.setValueOnly(offset, adaptivityValue);
2229 for (
auto it = signFlagsNode.cbeginValueOn(); it; ++it) {
2230 const Int16 flags = it.getValue();
2231 const unsigned char signs =
static_cast<unsigned char>(SIGNS &
int(flags));
2233 if ((flags & SEAM) || !sAdaptable[signs] || sEdgeGroupTable[signs][0] > 1) {
2235 mask.setActiveState(it.getCoord() & ~1u,
true);
2237 }
else if (flags & EDGES) {
2239 bool maskRegion =
false;
2241 ijk = it.getCoord();
2242 if (!pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2244 if (!maskRegion && flags & XEDGE) {
2246 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2248 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2250 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2254 if (!maskRegion && flags & YEDGE) {
2256 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2258 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2260 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2264 if (!maskRegion && flags & ZEDGE) {
2266 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2268 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2270 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2275 mask.setActiveState(it.getCoord() & ~1u,
true);
2282 for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
2283 for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
2284 for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
2285 if (!mask.isValueOn(ijk) && isNonManifold(inputAcc, ijk, mIsovalue, dim)) {
2286 mask.setActiveState(ijk,
true);
2297 gradientNode->setValuesOff();
2299 gradientNode.reset(
new Vec3sLeafNodeType());
2302 for (
auto it = signFlagsNode.cbeginValueOn(); it; ++it)
2304 ijk = it.getCoord();
2305 if (!mask.isValueOn(ijk & ~1u))
2310 if (invertGradientDir) {
2314 gradientNode->setValueOn(it.pos(), dir);
2321 for ( ; dim <= LeafDim; dim = dim << 1) {
2322 const unsigned coordMask = ~((dim << 1) - 1);
2323 for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
2324 for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
2325 for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
2327 adaptivity = adaptivityLeaf.getValue(ijk);
2329 if (mask.isValueOn(ijk)
2330 || isNonManifold(inputAcc, ijk, mIsovalue, dim)
2331 || (useGradients && !isMergeable(*gradientNode, ijk, dim, adaptivity)))
2333 mask.setActiveState(ijk & coordMask,
true);
2335 mergeVoxels(pointIndexNode, ijk, dim, regionId++);
2349 struct UniformPrimBuilder
2351 UniformPrimBuilder(): mIdx(0), mPolygonPool(nullptr) {}
2353 void init(
const size_t upperBound, PolygonPool& quadPool)
2355 mPolygonPool = &quadPool;
2356 mPolygonPool->resetQuads(upperBound);
2360 template<
typename IndexType>
2361 void addPrim(
const math::Vec4<IndexType>& verts,
bool reverse,
char flags = 0)
2364 mPolygonPool->quad(mIdx) = verts;
2366 Vec4I& quad = mPolygonPool->quad(mIdx);
2372 mPolygonPool->quadFlags(mIdx) =
flags;
2378 mPolygonPool->trimQuads(mIdx);
2383 PolygonPool* mPolygonPool;
2388 struct AdaptivePrimBuilder
2390 AdaptivePrimBuilder() : mQuadIdx(0), mTriangleIdx(0), mPolygonPool(nullptr) {}
2392 void init(
const size_t upperBound, PolygonPool& polygonPool)
2394 mPolygonPool = &polygonPool;
2395 mPolygonPool->resetQuads(upperBound);
2396 mPolygonPool->resetTriangles(upperBound);
2402 template<
typename IndexType>
2403 void addPrim(
const math::Vec4<IndexType>& verts,
bool reverse,
char flags = 0)
2405 if (verts[0] != verts[1] && verts[0] != verts[2] && verts[0] != verts[3]
2406 && verts[1] != verts[2] && verts[1] != verts[3] && verts[2] != verts[3]) {
2407 mPolygonPool->quadFlags(mQuadIdx) =
flags;
2408 addQuad(verts, reverse);
2410 verts[0] == verts[3] &&
2411 verts[1] != verts[2] &&
2412 verts[1] != verts[0] &&
2413 verts[2] != verts[0]) {
2414 mPolygonPool->triangleFlags(mTriangleIdx) =
flags;
2415 addTriangle(verts[0], verts[1], verts[2], reverse);
2417 verts[1] == verts[2] &&
2418 verts[0] != verts[3] &&
2419 verts[0] != verts[1] &&
2420 verts[3] != verts[1]) {
2421 mPolygonPool->triangleFlags(mTriangleIdx) =
flags;
2422 addTriangle(verts[0], verts[1], verts[3], reverse);
2424 verts[0] == verts[1] &&
2425 verts[2] != verts[3] &&
2426 verts[2] != verts[0] &&
2427 verts[3] != verts[0]) {
2428 mPolygonPool->triangleFlags(mTriangleIdx) =
flags;
2429 addTriangle(verts[0], verts[2], verts[3], reverse);
2431 verts[2] == verts[3] &&
2432 verts[0] != verts[1] &&
2433 verts[0] != verts[2] &&
2434 verts[1] != verts[2]) {
2435 mPolygonPool->triangleFlags(mTriangleIdx) =
flags;
2436 addTriangle(verts[0], verts[1], verts[2], reverse);
2443 mPolygonPool->trimQuads(mQuadIdx,
true);
2444 mPolygonPool->trimTrinagles(mTriangleIdx,
true);
2449 template<
typename IndexType>
2450 void addQuad(
const math::Vec4<IndexType>& verts,
bool reverse)
2453 mPolygonPool->quad(mQuadIdx) = verts;
2455 Vec4I& quad = mPolygonPool->quad(mQuadIdx);
2464 void addTriangle(
unsigned v0,
unsigned v1,
unsigned v2,
bool reverse)
2466 Vec3I& prim = mPolygonPool->triangle(mTriangleIdx);
2480 size_t mQuadIdx, mTriangleIdx;
2481 PolygonPool *mPolygonPool;
2485 template<
typename SignAccT,
typename IdxAccT,
typename PrimBuilder>
2488 bool invertSurfaceOrientation,
2493 const SignAccT& signAcc,
2494 const IdxAccT& idxAcc,
2495 PrimBuilder& mesher)
2500 const bool isActive = idxAcc.probeValue(ijk, v0);
2507 bool isInside = flags &
INSIDE;
2509 isInside = invertSurfaceOrientation ? !isInside : isInside;
2512 math::Vec4<IndexType> quad(0,0,0,0);
2514 if (flags & XEDGE) {
2516 quad[0] = v0 + offsets[0];
2520 bool activeValues = idxAcc.probeValue(coord, quad[1]);
2521 uint8_t cell = uint8_t(SIGNS & signAcc.getValue(coord));
2522 quad[1] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][5] - 1 : 0;
2526 activeValues = activeValues && idxAcc.probeValue(coord, quad[2]);
2527 cell = uint8_t(SIGNS & signAcc.getValue(coord));
2528 quad[2] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][7] - 1 : 0;
2532 activeValues = activeValues && idxAcc.probeValue(coord, quad[3]);
2533 cell = uint8_t(SIGNS & signAcc.getValue(coord));
2534 quad[3] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][3] - 1 : 0;
2537 mesher.addPrim(quad, isInside, tag[
bool(refFlags & XEDGE)]);
2544 if (flags & YEDGE) {
2546 quad[0] = v0 + offsets[1];
2550 bool activeValues = idxAcc.probeValue(coord, quad[1]);
2551 uint8_t cell = uint8_t(SIGNS & signAcc.getValue(coord));
2552 quad[1] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][12] - 1 : 0;
2556 activeValues = activeValues && idxAcc.probeValue(coord, quad[2]);
2557 cell = uint8_t(SIGNS & signAcc.getValue(coord));
2558 quad[2] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][11] - 1 : 0;
2562 activeValues = activeValues && idxAcc.probeValue(coord, quad[3]);
2563 cell = uint8_t(SIGNS & signAcc.getValue(coord));
2564 quad[3] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][10] - 1 : 0;
2567 mesher.addPrim(quad, isInside, tag[
bool(refFlags & YEDGE)]);
2574 if (flags & ZEDGE) {
2576 quad[0] = v0 + offsets[2];
2580 bool activeValues = idxAcc.probeValue(coord, quad[1]);
2581 uint8_t cell = uint8_t(SIGNS & signAcc.getValue(coord));
2582 quad[1] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][8] - 1 : 0;
2586 activeValues = activeValues && idxAcc.probeValue(coord, quad[2]);
2587 cell = uint8_t(SIGNS & signAcc.getValue(coord));
2588 quad[2] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][6] - 1 : 0;
2592 activeValues = activeValues && idxAcc.probeValue(coord, quad[3]);
2593 cell = uint8_t(SIGNS & signAcc.getValue(coord));
2594 quad[3] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][2] - 1 : 0;
2597 mesher.addPrim(quad, !isInside, tag[
bool(refFlags & ZEDGE)]);
2606 template<
typename InputTreeType>
2607 struct MaskTileBorders
2610 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
2613 MaskTileBorders(
const InputTreeType& inputTree, InputValueType iso,
2614 BoolTreeType& mask,
const Vec4i* tileArray)
2615 : mInputTree(&inputTree)
2619 , mTileArray(tileArray)
2623 MaskTileBorders(MaskTileBorders& rhs,
tbb::split)
2624 : mInputTree(rhs.mInputTree)
2625 , mIsovalue(rhs.mIsovalue)
2628 , mTileArray(rhs.mTileArray)
2632 void join(MaskTileBorders& rhs) {
mMask->merge(*rhs.mMask); }
2634 void operator()(
const tbb::blocked_range<size_t>&);
2637 InputTreeType
const *
const mInputTree;
2638 InputValueType
const mIsovalue;
2639 BoolTreeType mTempMask;
2640 BoolTreeType *
const mMask;
2641 Vec4i const *
const mTileArray;
2645 template<
typename InputTreeType>
2647 MaskTileBorders<InputTreeType>::operator()(
const tbb::blocked_range<size_t>& range)
2649 tree::ValueAccessor<const InputTreeType> inputTreeAcc(*mInputTree);
2654 for (
size_t n = range.begin(), N = range.end();
n !=
N; ++
n) {
2656 const Vec4i& tile = mTileArray[
n];
2658 bbox.min()[0] = tile[0];
2659 bbox.min()[1] = tile[1];
2660 bbox.min()[2] = tile[2];
2661 bbox.max() = bbox.min();
2662 bbox.max().offset(tile[3]);
2664 InputValueType
value = mInputTree->background();
2666 const bool isInside = isInsideValue(inputTreeAcc.getValue(bbox.min()), mIsovalue);
2667 const int valueDepth = inputTreeAcc.getValueDepth(bbox.min());
2675 bool processRegion =
true;
2676 if (valueDepth >= inputTreeAcc.getValueDepth(nijk)) {
2677 processRegion = isInside != isInsideValue(inputTreeAcc.getValue(nijk), mIsovalue);
2680 if (processRegion) {
2683 region.min()[0] = region.max()[0] = ijk[0];
2684 mMask->fill(region,
false);
2691 processRegion =
true;
2692 if (valueDepth >= inputTreeAcc.getValueDepth(ijk)) {
2693 processRegion = (!inputTreeAcc.probeValue(ijk, value)
2694 && isInside != isInsideValue(value, mIsovalue));
2697 if (processRegion) {
2700 region.min()[0] = region.max()[0] = ijk[0];
2701 mMask->fill(region,
false);
2711 processRegion =
true;
2712 if (valueDepth >= inputTreeAcc.getValueDepth(nijk)) {
2713 processRegion = isInside != isInsideValue(inputTreeAcc.getValue(nijk), mIsovalue);
2716 if (processRegion) {
2719 region.min()[1] = region.max()[1] = ijk[1];
2720 mMask->fill(region,
false);
2727 processRegion =
true;
2728 if (valueDepth >= inputTreeAcc.getValueDepth(ijk)) {
2729 processRegion = (!inputTreeAcc.probeValue(ijk, value)
2730 && isInside != isInsideValue(value, mIsovalue));
2733 if (processRegion) {
2736 region.min()[1] = region.max()[1] = ijk[1];
2737 mMask->fill(region,
false);
2747 processRegion =
true;
2748 if (valueDepth >= inputTreeAcc.getValueDepth(nijk)) {
2749 processRegion = isInside != isInsideValue(inputTreeAcc.getValue(nijk), mIsovalue);
2752 if (processRegion) {
2755 region.min()[2] = region.max()[2] = ijk[2];
2756 mMask->fill(region,
false);
2762 processRegion =
true;
2763 if (valueDepth >= inputTreeAcc.getValueDepth(ijk)) {
2764 processRegion = (!inputTreeAcc.probeValue(ijk, value)
2765 && isInside != isInsideValue(value, mIsovalue));
2768 if (processRegion) {
2771 region.min()[2] = region.max()[2] = ijk[2];
2772 mMask->fill(region,
false);
2778 template<
typename InputTreeType>
2780 maskActiveTileBorders(
const InputTreeType& inputTree,
2782 typename InputTreeType::template ValueConverter<bool>::Type& mask)
2784 typename InputTreeType::ValueOnCIter tileIter(inputTree);
2785 tileIter.setMaxDepth(InputTreeType::ValueOnCIter::LEAF_DEPTH - 1);
2787 size_t tileCount = 0;
2788 for ( ; tileIter; ++tileIter) {
2792 if (tileCount > 0) {
2793 std::unique_ptr<Vec4i[]> tiles(
new Vec4i[tileCount]);
2798 tileIter = inputTree.cbeginValueOn();
2799 tileIter.setMaxDepth(InputTreeType::ValueOnCIter::LEAF_DEPTH - 1);
2801 for (; tileIter; ++tileIter) {
2802 Vec4i& tile = tiles[index++];
2803 tileIter.getBoundingBox(bbox);
2804 tile[0] = bbox.min()[0];
2805 tile[1] = bbox.min()[1];
2806 tile[2] = bbox.min()[2];
2807 tile[3] = bbox.max()[0] - bbox.min()[0];
2810 MaskTileBorders<InputTreeType> op(inputTree, iso, mask, tiles.get());
2811 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, tileCount), op);
2823 PointListCopy(
const PointList& pointsIn, std::vector<Vec3s>& pointsOut)
2824 : mPointsIn(pointsIn) , mPointsOut(pointsOut)
2828 void operator()(
const tbb::blocked_range<size_t>& range)
const
2830 for (
size_t n = range.begin();
n < range.end(); ++
n) {
2831 mPointsOut[
n] = mPointsIn[
n];
2837 std::vector<Vec3s>& mPointsOut;
2845 struct LeafNodeVoxelOffsets
2847 using IndexVector = std::vector<Index>;
2849 template<
typename LeafNodeType>
2850 void constructOffsetList();
2853 const IndexVector& core()
const {
return mCore; }
2857 const IndexVector& minX()
const {
return mMinX; }
2860 const IndexVector& maxX()
const {
return mMaxX; }
2864 const IndexVector& minY()
const {
return mMinY; }
2867 const IndexVector& maxY()
const {
return mMaxY; }
2871 const IndexVector& minZ()
const {
return mMinZ; }
2874 const IndexVector& maxZ()
const {
return mMaxZ; }
2878 const IndexVector& internalNeighborsX()
const {
return mInternalNeighborsX; }
2881 const IndexVector& internalNeighborsY()
const {
return mInternalNeighborsY; }
2884 const IndexVector& internalNeighborsZ()
const {
return mInternalNeighborsZ; }
2888 IndexVector mCore, mMinX, mMaxX, mMinY, mMaxY, mMinZ, mMaxZ,
2889 mInternalNeighborsX, mInternalNeighborsY, mInternalNeighborsZ;
2893 template<
typename LeafNodeType>
2895 LeafNodeVoxelOffsets::constructOffsetList()
2899 mCore.reserve((LeafNodeType::DIM - 2) * (LeafNodeType::DIM - 2));
2901 for (
Index x = 1;
x < (LeafNodeType::DIM - 1); ++
x) {
2902 const Index offsetX =
x << (2 * LeafNodeType::LOG2DIM);
2903 for (
Index y = 1;
y < (LeafNodeType::DIM - 1); ++
y) {
2904 const Index offsetXY = offsetX + (
y << LeafNodeType::LOG2DIM);
2905 for (
Index z = 1;
z < (LeafNodeType::DIM - 1); ++
z) {
2906 mCore.push_back(offsetXY +
z);
2912 mInternalNeighborsX.clear();
2913 mInternalNeighborsX.reserve(
LeafNodeType::SIZE - (LeafNodeType::DIM * LeafNodeType::DIM));
2915 for (
Index x = 0;
x < (LeafNodeType::DIM - 1); ++
x) {
2916 const Index offsetX =
x << (2 * LeafNodeType::LOG2DIM);
2917 for (
Index y = 0;
y < LeafNodeType::DIM; ++
y) {
2918 const Index offsetXY = offsetX + (
y << LeafNodeType::LOG2DIM);
2919 for (
Index z = 0;
z < LeafNodeType::DIM; ++
z) {
2920 mInternalNeighborsX.push_back(offsetXY +
z);
2926 mInternalNeighborsY.clear();
2927 mInternalNeighborsY.reserve(
LeafNodeType::SIZE - (LeafNodeType::DIM * LeafNodeType::DIM));
2929 for (
Index x = 0;
x < LeafNodeType::DIM; ++
x) {
2930 const Index offsetX =
x << (2 * LeafNodeType::LOG2DIM);
2931 for (
Index y = 0;
y < (LeafNodeType::DIM - 1); ++
y) {
2932 const Index offsetXY = offsetX + (
y << LeafNodeType::LOG2DIM);
2933 for (
Index z = 0;
z < LeafNodeType::DIM; ++
z) {
2934 mInternalNeighborsY.push_back(offsetXY +
z);
2940 mInternalNeighborsZ.clear();
2941 mInternalNeighborsZ.reserve(
LeafNodeType::SIZE - (LeafNodeType::DIM * LeafNodeType::DIM));
2943 for (
Index x = 0;
x < LeafNodeType::DIM; ++
x) {
2944 const Index offsetX =
x << (2 * LeafNodeType::LOG2DIM);
2945 for (
Index y = 0;
y < LeafNodeType::DIM; ++
y) {
2946 const Index offsetXY = offsetX + (
y << LeafNodeType::LOG2DIM);
2947 for (
Index z = 0;
z < (LeafNodeType::DIM - 1); ++
z) {
2948 mInternalNeighborsZ.push_back(offsetXY +
z);
2955 mMinX.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
2957 for (
Index y = 0;
y < LeafNodeType::DIM; ++
y) {
2958 const Index offsetXY = (
y << LeafNodeType::LOG2DIM);
2959 for (
Index z = 0;
z < LeafNodeType::DIM; ++
z) {
2960 mMinX.push_back(offsetXY +
z);
2967 mMaxX.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
2969 const Index offsetX = (LeafNodeType::DIM - 1) << (2 * LeafNodeType::LOG2DIM);
2970 for (
Index y = 0;
y < LeafNodeType::DIM; ++
y) {
2971 const Index offsetXY = offsetX + (
y << LeafNodeType::LOG2DIM);
2972 for (
Index z = 0;
z < LeafNodeType::DIM; ++
z) {
2973 mMaxX.push_back(offsetXY +
z);
2980 mMinY.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
2982 for (
Index x = 0;
x < LeafNodeType::DIM; ++
x) {
2983 const Index offsetX =
x << (2 * LeafNodeType::LOG2DIM);
2984 for (
Index z = 0;
z < (LeafNodeType::DIM - 1); ++
z) {
2985 mMinY.push_back(offsetX +
z);
2992 mMaxY.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
2994 const Index offsetY = (LeafNodeType::DIM - 1) << LeafNodeType::LOG2DIM;
2995 for (
Index x = 0;
x < LeafNodeType::DIM; ++
x) {
2996 const Index offsetX =
x << (2 * LeafNodeType::LOG2DIM);
2997 for (
Index z = 0;
z < (LeafNodeType::DIM - 1); ++
z) {
2998 mMaxY.push_back(offsetX + offsetY +
z);
3005 mMinZ.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
3007 for (
Index x = 0;
x < LeafNodeType::DIM; ++
x) {
3008 const Index offsetX =
x << (2 * LeafNodeType::LOG2DIM);
3009 for (
Index y = 0;
y < LeafNodeType::DIM; ++
y) {
3010 const Index offsetXY = offsetX + (
y << LeafNodeType::LOG2DIM);
3011 mMinZ.push_back(offsetXY);
3018 mMaxZ.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
3020 for (
Index x = 0;
x < LeafNodeType::DIM; ++
x) {
3021 const Index offsetX =
x << (2 * LeafNodeType::LOG2DIM);
3022 for (
Index y = 0;
y < LeafNodeType::DIM; ++
y) {
3023 const Index offsetXY = offsetX + (
y << LeafNodeType::LOG2DIM);
3024 mMaxZ.push_back(offsetXY + (LeafNodeType::DIM - 1));
3035 template<
typename AccessorT,
int _AXIS>
3036 struct VoxelEdgeAccessor
3038 enum {
AXIS = _AXIS };
3041 VoxelEdgeAccessor(AccessorT& _acc) : acc(_acc) {}
3043 void set(Coord ijk) {
3045 acc.setActiveState(ijk);
3047 acc.setActiveState(ijk);
3049 acc.setActiveState(ijk);
3051 acc.setActiveState(ijk);
3052 }
else if (_AXIS == 1) {
3053 acc.setActiveState(ijk);
3055 acc.setActiveState(ijk);
3057 acc.setActiveState(ijk);
3059 acc.setActiveState(ijk);
3061 acc.setActiveState(ijk);
3063 acc.setActiveState(ijk);
3065 acc.setActiveState(ijk);
3067 acc.setActiveState(ijk);
3076 template<
typename VoxelEdgeAcc,
typename LeafNodeT>
3078 evalInternalVoxelEdges(VoxelEdgeAcc& edgeAcc,
3079 const LeafNodeT& leafnode,
3080 const LeafNodeVoxelOffsets& voxels,
3084 const std::vector<Index>* offsets = &voxels.internalNeighborsZ();
3086 if (VoxelEdgeAcc::AXIS == 0) {
3087 nvo = LeafNodeT::DIM * LeafNodeT::DIM;
3088 offsets = &voxels.internalNeighborsX();
3089 }
else if (VoxelEdgeAcc::AXIS == 1) {
3090 nvo = LeafNodeT::DIM;
3091 offsets = &voxels.internalNeighborsY();
3094 const LeafBufferAccessor<LeafNodeT> lhsAcc(leafnode);
3096 for (
size_t n = 0, N = offsets->size();
n <
N; ++
n) {
3097 const Index& pos = (*offsets)[
n];
3098 const bool isActive = leafnode.isValueOn(pos) || leafnode.isValueOn(pos + nvo);
3099 if (isActive && (isInsideValue(lhsAcc.get(pos), iso) !=
3100 isInsideValue(lhsAcc.get((pos + nvo)), iso))) {
3101 edgeAcc.set(leafnode.offsetToGlobalCoord(pos));
3110 template<
typename LeafNodeT,
typename TreeAcc,
typename VoxelEdgeAcc>
3112 evalExternalVoxelEdges(VoxelEdgeAcc& edgeAcc,
3114 const LeafNodeT& lhsNode,
3115 const LeafNodeVoxelOffsets& voxels,
3118 const std::vector<Index>* lhsOffsets = &voxels.maxX();
3119 const std::vector<Index>* rhsOffsets = &voxels.minX();
3120 Coord ijk = lhsNode.origin();
3122 if (VoxelEdgeAcc::AXIS == 0) {
3123 ijk[0] += LeafNodeT::DIM;
3124 }
else if (VoxelEdgeAcc::AXIS == 1) {
3125 ijk[1] += LeafNodeT::DIM;
3126 lhsOffsets = &voxels.maxY();
3127 rhsOffsets = &voxels.minY();
3128 }
else if (VoxelEdgeAcc::AXIS == 2) {
3129 ijk[2] += LeafNodeT::DIM;
3130 lhsOffsets = &voxels.maxZ();
3131 rhsOffsets = &voxels.minZ();
3135 const LeafNodeT* rhsNodePt = acc.probeConstLeaf(ijk);
3137 const LeafBufferAccessor<LeafNodeT> lhsAcc(lhsNode);
3140 const LeafBufferAccessor<LeafNodeT> rhsAcc(*rhsNodePt);
3142 for (
size_t n = 0, N = lhsOffsets->size();
n <
N; ++
n) {
3143 const Index& pos = (*lhsOffsets)[
n];
3144 bool isActive = lhsNode.isValueOn(pos) || rhsNodePt->isValueOn((*rhsOffsets)[
n]);
3145 if (isActive && (isInsideValue(lhsAcc.get(pos), iso) !=
3146 isInsideValue(rhsAcc.get((*rhsOffsets)[n]), iso))) {
3147 edgeAcc.set(lhsNode.offsetToGlobalCoord(pos));
3150 }
else if (!acc.probeValue(ijk, value)) {
3151 const bool inside = isInsideValue(value, iso);
3152 for (
size_t n = 0, N = lhsOffsets->size(); n <
N; ++
n) {
3153 const Index& pos = (*lhsOffsets)[
n];
3154 if (lhsNode.isValueOn(pos) && (inside != isInsideValue(lhsAcc.get(pos), iso))) {
3155 edgeAcc.set(lhsNode.offsetToGlobalCoord(pos));
3165 template<
typename LeafNodeT,
typename TreeAcc,
typename VoxelEdgeAcc>
3167 evalExternalVoxelEdgesInv(VoxelEdgeAcc& edgeAcc,
3169 const LeafNodeT& leafnode,
3170 const LeafNodeVoxelOffsets& voxels,
3173 Coord ijk = leafnode.origin();
3174 if (VoxelEdgeAcc::AXIS == 0) --ijk[0];
3175 else if (VoxelEdgeAcc::AXIS == 1) --ijk[1];
3176 else if (VoxelEdgeAcc::AXIS == 2) --ijk[2];
3179 if (!acc.probeConstLeaf(ijk) && !acc.probeValue(ijk, value)) {
3181 const std::vector<Index>* offsets = &voxels.internalNeighborsX();
3182 if (VoxelEdgeAcc::AXIS == 1) offsets = &voxels.internalNeighborsY();
3183 else if (VoxelEdgeAcc::AXIS == 2) offsets = &voxels.internalNeighborsZ();
3185 const LeafBufferAccessor<LeafNodeT> lhsAcc(leafnode);
3187 const bool inside = isInsideValue(value, iso);
3188 for (
size_t n = 0, N = offsets->size(); n <
N; ++
n) {
3190 const Index& pos = (*offsets)[
n];
3191 if (leafnode.isValueOn(pos)
3192 && (inside != isInsideValue(lhsAcc.get(pos), iso)))
3194 ijk = leafnode.offsetToGlobalCoord(pos);
3195 if (VoxelEdgeAcc::AXIS == 0) --ijk[0];
3196 else if (VoxelEdgeAcc::AXIS == 1) --ijk[1];
3197 else if (VoxelEdgeAcc::AXIS == 2) --ijk[2];
3206 template<
typename InputTreeType>
3207 struct IdentifyIntersectingVoxels
3209 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
3212 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
3214 IdentifyIntersectingVoxels(
3215 const InputTreeType& inputTree,
3216 const std::vector<const InputLeafNodeType*>& inputLeafNodes,
3217 const LeafNodeVoxelOffsets& offsets,
3218 BoolTreeType& intersectionTree,
3219 InputValueType iso);
3221 IdentifyIntersectingVoxels(IdentifyIntersectingVoxels&,
tbb::split);
3222 void operator()(
const tbb::blocked_range<size_t>&);
3223 void join(
const IdentifyIntersectingVoxels& rhs) {
3224 mIntersectionAccessor.tree().merge(rhs.mIntersectionAccessor.tree());
3228 tree::ValueAccessor<const InputTreeType> mInputAccessor;
3229 InputLeafNodeType
const *
const *
const mInputNodes;
3231 BoolTreeType mIntersectionTree;
3232 tree::ValueAccessor<BoolTreeType> mIntersectionAccessor;
3234 const LeafNodeVoxelOffsets& mOffsets;
3235 const InputValueType mIsovalue;
3239 template<
typename InputTreeType>
3240 IdentifyIntersectingVoxels<InputTreeType>::IdentifyIntersectingVoxels(
3241 const InputTreeType& inputTree,
3242 const std::vector<const InputLeafNodeType*>& inputLeafNodes,
3243 const LeafNodeVoxelOffsets& offsets,
3244 BoolTreeType& intersectionTree,
3246 : mInputAccessor(inputTree)
3247 , mInputNodes(inputLeafNodes.data())
3248 , mIntersectionTree(false)
3249 , mIntersectionAccessor(intersectionTree)
3256 template<
typename InputTreeType>
3257 IdentifyIntersectingVoxels<InputTreeType>::IdentifyIntersectingVoxels(
3259 : mInputAccessor(rhs.mInputAccessor.tree())
3260 , mInputNodes(rhs.mInputNodes)
3261 , mIntersectionTree(false)
3262 , mIntersectionAccessor(mIntersectionTree)
3263 , mOffsets(rhs.mOffsets)
3264 , mIsovalue(rhs.mIsovalue)
3269 template<
typename InputTreeType>
3271 IdentifyIntersectingVoxels<InputTreeType>::operator()(
const tbb::blocked_range<size_t>& range)
3273 VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 0> xEdgeAcc(mIntersectionAccessor);
3274 VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 1> yEdgeAcc(mIntersectionAccessor);
3275 VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 2> zEdgeAcc(mIntersectionAccessor);
3277 for (
size_t n = range.begin(); n != range.end(); ++
n) {
3279 const InputLeafNodeType& node = *mInputNodes[
n];
3282 evalInternalVoxelEdges(xEdgeAcc, node, mOffsets, mIsovalue);
3284 evalInternalVoxelEdges(yEdgeAcc, node, mOffsets, mIsovalue);
3286 evalInternalVoxelEdges(zEdgeAcc, node, mOffsets, mIsovalue);
3289 evalExternalVoxelEdges(xEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3291 evalExternalVoxelEdges(yEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3293 evalExternalVoxelEdges(zEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3299 evalExternalVoxelEdgesInv(xEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3301 evalExternalVoxelEdgesInv(yEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3303 evalExternalVoxelEdgesInv(zEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3308 template<
typename InputTreeType>
3310 identifySurfaceIntersectingVoxels(
3311 typename InputTreeType::template ValueConverter<bool>::Type& intersectionTree,
3312 const InputTreeType& inputTree,
3315 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
3317 std::vector<const InputLeafNodeType*> inputLeafNodes;
3318 inputTree.getNodes(inputLeafNodes);
3321 offsets.constructOffsetList<InputLeafNodeType>();
3323 IdentifyIntersectingVoxels<InputTreeType> op(
3324 inputTree, inputLeafNodes, offsets, intersectionTree, isovalue);
3326 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, inputLeafNodes.size()), op);
3328 maskActiveTileBorders(inputTree, isovalue, intersectionTree);
3335 template<
typename InputTreeType>
3336 struct MaskIntersectingVoxels
3338 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
3341 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
3342 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
3344 MaskIntersectingVoxels(
3345 const InputTreeType& inputTree,
3346 const std::vector<BoolLeafNodeType*>& nodes,
3347 BoolTreeType& intersectionTree,
3348 InputValueType iso);
3350 MaskIntersectingVoxels(MaskIntersectingVoxels&,
tbb::split);
3351 void operator()(
const tbb::blocked_range<size_t>&);
3352 void join(
const MaskIntersectingVoxels& rhs) {
3353 mIntersectionAccessor.tree().merge(rhs.mIntersectionAccessor.tree());
3357 tree::ValueAccessor<const InputTreeType> mInputAccessor;
3358 BoolLeafNodeType
const *
const *
const mNodes;
3360 BoolTreeType mIntersectionTree;
3361 tree::ValueAccessor<BoolTreeType> mIntersectionAccessor;
3363 const InputValueType mIsovalue;
3367 template<
typename InputTreeType>
3368 MaskIntersectingVoxels<InputTreeType>::MaskIntersectingVoxels(
3369 const InputTreeType& inputTree,
3370 const std::vector<BoolLeafNodeType*>& nodes,
3371 BoolTreeType& intersectionTree,
3373 : mInputAccessor(inputTree)
3374 , mNodes(nodes.data())
3375 , mIntersectionTree(false)
3376 , mIntersectionAccessor(intersectionTree)
3382 template<
typename InputTreeType>
3383 MaskIntersectingVoxels<InputTreeType>::MaskIntersectingVoxels(
3385 : mInputAccessor(rhs.mInputAccessor.tree())
3386 , mNodes(rhs.mNodes)
3387 , mIntersectionTree(false)
3388 , mIntersectionAccessor(mIntersectionTree)
3389 , mIsovalue(rhs.mIsovalue)
3394 template<
typename InputTreeType>
3396 MaskIntersectingVoxels<InputTreeType>::operator()(
const tbb::blocked_range<size_t>& range)
3398 VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 0> xEdgeAcc(mIntersectionAccessor);
3399 VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 1> yEdgeAcc(mIntersectionAccessor);
3400 VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 2> zEdgeAcc(mIntersectionAccessor);
3404 for (
size_t n = range.begin(); n != range.end(); ++
n) {
3406 const BoolLeafNodeType& node = *mNodes[
n];
3408 for (
auto it = node.cbeginValueOn(); it; ++it) {
3410 if (!it.getValue()) {
3412 ijk = it.getCoord();
3414 const bool inside = isInsideValue(mInputAccessor.getValue(ijk), mIsovalue);
3416 if (inside != isInsideValue(mInputAccessor.getValue(ijk.offsetBy(1, 0, 0)), mIsovalue)) {
3420 if (inside != isInsideValue(mInputAccessor.getValue(ijk.offsetBy(0, 1, 0)), mIsovalue)) {
3424 if (inside != isInsideValue(mInputAccessor.getValue(ijk.offsetBy(0, 0, 1)), mIsovalue)) {
3433 template<
typename BoolTreeType>
3434 struct MaskBorderVoxels
3436 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
3438 MaskBorderVoxels(
const BoolTreeType& maskTree,
3439 const std::vector<BoolLeafNodeType*>& maskNodes,
3440 BoolTreeType& borderTree)
3441 : mMaskTree(&maskTree)
3442 , mMaskNodes(maskNodes.data())
3443 , mTmpBorderTree(false)
3444 , mBorderTree(&borderTree) {}
3446 MaskBorderVoxels(MaskBorderVoxels& rhs,
tbb::split)
3447 : mMaskTree(rhs.mMaskTree)
3448 , mMaskNodes(rhs.mMaskNodes)
3449 , mTmpBorderTree(false)
3450 , mBorderTree(&mTmpBorderTree) {}
3452 void join(MaskBorderVoxels& rhs) { mBorderTree->merge(*rhs.mBorderTree); }
3454 void operator()(
const tbb::blocked_range<size_t>& range)
3456 tree::ValueAccessor<const BoolTreeType> maskAcc(*mMaskTree);
3457 tree::ValueAccessor<BoolTreeType> borderAcc(*mBorderTree);
3460 for (
size_t n = range.begin(); n != range.end(); ++
n) {
3462 const BoolLeafNodeType& node = *mMaskNodes[
n];
3464 for (
auto it = node.cbeginValueOn(); it; ++it) {
3466 ijk = it.getCoord();
3468 const bool lhs = it.getValue();
3471 bool isEdgeVoxel =
false;
3474 isEdgeVoxel = (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3477 isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3480 isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3483 isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3487 isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3490 isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3493 isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3497 borderAcc.setValue(ijk,
true);
3504 BoolTreeType
const *
const mMaskTree;
3505 BoolLeafNodeType
const *
const *
const mMaskNodes;
3507 BoolTreeType mTmpBorderTree;
3508 BoolTreeType *
const mBorderTree;
3512 template<
typename BoolTreeType>
3513 struct SyncMaskValues
3515 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
3517 SyncMaskValues(
const std::vector<BoolLeafNodeType*>& nodes,
const BoolTreeType& mask)
3518 : mNodes(nodes.data())
3519 , mMaskTree(&mask) {}
3521 void operator()(
const tbb::blocked_range<size_t>& range)
const
3523 using ValueOnIter =
typename BoolLeafNodeType::ValueOnIter;
3525 tree::ValueAccessor<const BoolTreeType> maskTreeAcc(*mMaskTree);
3527 for (
size_t n = range.begin(), N = range.end(); n !=
N; ++
n) {
3529 BoolLeafNodeType& node = *mNodes[
n];
3531 const BoolLeafNodeType * maskNode = maskTreeAcc.probeConstLeaf(node.origin());
3534 for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3535 const Index pos = it.pos();
3536 if (maskNode->getValue(pos)) {
3537 node.setValueOnly(pos,
true);
3545 BoolLeafNodeType *
const *
const mNodes;
3546 BoolTreeType
const *
const mMaskTree;
3553 template<
typename BoolTreeType>
3556 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
3558 MaskSurface(
const std::vector<BoolLeafNodeType*>& nodes,
3559 const BoolTreeType& mask,
3560 const math::Transform& inputTransform,
3561 const math::Transform& maskTransform,
3563 : mNodes(nodes.data())
3565 , mInputTransform(inputTransform)
3566 , mMaskTransform(maskTransform)
3567 , mMatchingTransforms(mInputTransform == mMaskTransform)
3568 , mInvertMask(invert) {}
3570 void operator()(
const tbb::blocked_range<size_t>& range)
const
3572 using ValueOnIter =
typename BoolLeafNodeType::ValueOnIter;
3574 tree::ValueAccessor<const BoolTreeType> maskTreeAcc(*mMaskTree);
3576 for (
size_t n = range.begin(), N = range.end(); n !=
N; ++
n) {
3578 BoolLeafNodeType& node = *mNodes[
n];
3580 if (mMatchingTransforms) {
3582 const BoolLeafNodeType * maskNode = maskTreeAcc.probeConstLeaf(node.origin());
3586 for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3587 const Index pos = it.pos();
3588 if (maskNode->isValueOn(pos) == mInvertMask) {
3589 node.setValueOnly(pos,
true);
3595 if (maskTreeAcc.isValueOn(node.origin()) == mInvertMask) {
3596 for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3597 node.setValueOnly(it.pos(),
true);
3607 for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3609 ijk = mMaskTransform.worldToIndexCellCentered(
3610 mInputTransform.indexToWorld(it.getCoord()));
3612 if (maskTreeAcc.isValueOn(ijk) == mInvertMask) {
3613 node.setValueOnly(it.pos(),
true);
3622 BoolLeafNodeType *
const *
const mNodes;
3623 BoolTreeType
const *
const mMaskTree;
3624 const math::Transform& mInputTransform;
3625 const math::Transform& mMaskTransform;
3626 const bool mMatchingTransforms;
3627 const bool mInvertMask;
3631 template<
typename InputGr
idType>
3634 typename InputGridType::TreeType::template ValueConverter<bool>::Type& intersectionTree,
3635 typename InputGridType::TreeType::template ValueConverter<bool>::Type& borderTree,
3636 const InputGridType& inputGrid,
3638 const bool invertMask,
3641 using InputTreeType =
typename InputGridType::TreeType;
3642 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
3643 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
3644 using BoolGridType = Grid<BoolTreeType>;
3646 if (!maskGrid)
return;
3647 if (maskGrid->type() != BoolGridType::gridType())
return;
3649 const math::Transform&
transform = inputGrid.transform();
3650 const InputTreeType& inputTree = inputGrid.tree();
3652 const BoolGridType * surfaceMask =
static_cast<const BoolGridType*
>(maskGrid.get());
3654 const BoolTreeType& maskTree = surfaceMask->tree();
3655 const math::Transform& maskTransform = surfaceMask->transform();
3659 std::vector<BoolLeafNodeType*> intersectionLeafNodes;
3660 intersectionTree.getNodes(intersectionLeafNodes);
3662 const tbb::blocked_range<size_t> intersectionRange(0, intersectionLeafNodes.size());
3665 MaskSurface<BoolTreeType>(
3666 intersectionLeafNodes, maskTree, transform, maskTransform, invertMask));
3671 MaskBorderVoxels<BoolTreeType> borderOp(
3672 intersectionTree, intersectionLeafNodes, borderTree);
3673 tbb::parallel_reduce(intersectionRange, borderOp);
3678 BoolTreeType tmpIntersectionTree(
false);
3680 MaskIntersectingVoxels<InputTreeType> op(
3681 inputTree, intersectionLeafNodes, tmpIntersectionTree, isovalue);
3683 tbb::parallel_reduce(intersectionRange, op);
3685 std::vector<BoolLeafNodeType*> tmpIntersectionLeafNodes;
3686 tmpIntersectionTree.getNodes(tmpIntersectionLeafNodes);
3688 tbb::parallel_for(tbb::blocked_range<size_t>(0, tmpIntersectionLeafNodes.size()),
3689 SyncMaskValues<BoolTreeType>(tmpIntersectionLeafNodes, intersectionTree));
3691 intersectionTree.clear();
3692 intersectionTree.merge(tmpIntersectionTree);
3699 template<
typename InputTreeType>
3700 struct ComputeAuxiliaryData
3702 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
3705 using BoolLeafNodeType = tree::LeafNode<bool, InputLeafNodeType::LOG2DIM>;
3711 ComputeAuxiliaryData(
const InputTreeType& inputTree,
3712 const std::vector<const BoolLeafNodeType*>& intersectionLeafNodes,
3713 Int16TreeType& signFlagsTree,
3714 Index32TreeType& pointIndexTree,
3715 InputValueType iso);
3717 ComputeAuxiliaryData(ComputeAuxiliaryData&,
tbb::split);
3718 void operator()(
const tbb::blocked_range<size_t>&);
3719 void join(
const ComputeAuxiliaryData& rhs) {
3720 mSignFlagsAccessor.tree().merge(rhs.mSignFlagsAccessor.tree());
3721 mPointIndexAccessor.tree().merge(rhs.mPointIndexAccessor.tree());
3725 tree::ValueAccessor<const InputTreeType> mInputAccessor;
3726 BoolLeafNodeType
const *
const *
const mIntersectionNodes;
3728 Int16TreeType mSignFlagsTree;
3729 tree::ValueAccessor<Int16TreeType> mSignFlagsAccessor;
3730 Index32TreeType mPointIndexTree;
3731 tree::ValueAccessor<Index32TreeType> mPointIndexAccessor;
3733 const InputValueType mIsovalue;
3737 template<
typename InputTreeType>
3738 ComputeAuxiliaryData<InputTreeType>::ComputeAuxiliaryData(
3739 const InputTreeType& inputTree,
3740 const std::vector<const BoolLeafNodeType*>& intersectionLeafNodes,
3741 Int16TreeType& signFlagsTree,
3742 Index32TreeType& pointIndexTree,
3744 : mInputAccessor(inputTree)
3745 , mIntersectionNodes(intersectionLeafNodes.data())
3747 , mSignFlagsAccessor(signFlagsTree)
3748 , mPointIndexTree(std::numeric_limits<
Index32>::
max())
3749 , mPointIndexAccessor(pointIndexTree)
3756 template<
typename InputTreeType>
3757 ComputeAuxiliaryData<InputTreeType>::ComputeAuxiliaryData(ComputeAuxiliaryData& rhs,
tbb::split)
3758 : mInputAccessor(rhs.mInputAccessor.tree())
3759 , mIntersectionNodes(rhs.mIntersectionNodes)
3761 , mSignFlagsAccessor(mSignFlagsTree)
3762 , mPointIndexTree(std::numeric_limits<
Index32>::
max())
3763 , mPointIndexAccessor(mPointIndexTree)
3764 , mIsovalue(rhs.mIsovalue)
3769 template<
typename InputTreeType>
3771 ComputeAuxiliaryData<InputTreeType>::operator()(
const tbb::blocked_range<size_t>& range)
3773 using Int16LeafNodeType =
typename Int16TreeType::LeafNodeType;
3776 std::array<InputValueType, 8> cellVertexValues;
3777 std::unique_ptr<Int16LeafNodeType> signsNodePt(
nullptr);
3779 for (
size_t n = range.begin(), N = range.end(); n !=
N; ++
n) {
3781 const BoolLeafNodeType& maskNode = *mIntersectionNodes[
n];
3782 const Coord& origin = maskNode.origin();
3784 const InputLeafNodeType* leafPt = mInputAccessor.probeConstLeaf(origin);
3786 if (!signsNodePt.get()) signsNodePt.reset(
new Int16LeafNodeType(origin, 0));
3787 else signsNodePt->setOrigin(origin);
3789 bool updatedNode =
false;
3791 for (
auto it = maskNode.cbeginValueOn(); it; ++it) {
3793 const Index pos = it.pos();
3794 ijk = BoolLeafNodeType::offsetToLocalCoord(pos);
3796 const bool inclusiveCell = leafPt && isInternalLeafCoord<InputLeafNodeType>(ijk);
3798 if (inclusiveCell) getCellVertexValues(*leafPt, pos, cellVertexValues);
3799 else getCellVertexValues(mInputAccessor, origin + ijk, cellVertexValues);
3801 uint8_t signFlags = computeSignFlags(cellVertexValues, mIsovalue);
3803 if (signFlags != 0 && signFlags != 0xFF) {
3805 const bool inside = signFlags & 0x1;
3807 int edgeFlags = inside ? INSIDE : 0;
3809 if (!it.getValue()) {
3810 edgeFlags |= inside != ((signFlags & 0x02) != 0) ? XEDGE : 0;
3811 edgeFlags |= inside != ((signFlags & 0x10) != 0) ? YEDGE : 0;
3812 edgeFlags |= inside != ((signFlags & 0x08) != 0) ? ZEDGE : 0;
3815 const uint8_t ambiguousCellFlags = sAmbiguousFace[signFlags];
3816 if (ambiguousCellFlags != 0) {
3817 correctCellSigns(signFlags, ambiguousCellFlags, mInputAccessor,
3818 origin + ijk, mIsovalue);
3821 edgeFlags |=
int(signFlags);
3823 signsNodePt->setValueOn(pos,
Int16(edgeFlags));
3829 typename Index32TreeType::LeafNodeType* idxNode =
3830 mPointIndexAccessor.touchLeaf(origin);
3831 idxNode->topologyUnion(*signsNodePt);
3834 auto*
const idxData = idxNode->buffer().data();
3835 for (
auto it = idxNode->beginValueOn(); it; ++it) {
3836 idxData[it.pos()] = 0;
3839 mSignFlagsAccessor.addLeaf(signsNodePt.release());
3845 template<
typename InputTreeType>
3847 computeAuxiliaryData(
3850 const typename InputTreeType::template ValueConverter<bool>::Type& intersectionTree,
3851 const InputTreeType& inputTree,
3854 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
3855 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
3857 std::vector<const BoolLeafNodeType*> intersectionLeafNodes;
3858 intersectionTree.getNodes(intersectionLeafNodes);
3860 ComputeAuxiliaryData<InputTreeType> op(
3861 inputTree, intersectionLeafNodes, signFlagsTree, pointIndexTree, isovalue);
3863 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, intersectionLeafNodes.size()), op);
3870 template<Index32 LeafNodeLog2Dim>
3871 struct LeafNodePointCount
3873 using Int16LeafNodeType = tree::LeafNode<Int16, LeafNodeLog2Dim>;
3875 LeafNodePointCount(
const std::vector<Int16LeafNodeType*>& leafNodes,
3876 std::unique_ptr<
Index32[]>& leafNodeCount)
3877 : mLeafNodes(leafNodes.data())
3878 , mData(leafNodeCount.
get()) {}
3880 void operator()(
const tbb::blocked_range<size_t>& range)
const {
3882 for (
size_t n = range.begin(), N = range.end(); n !=
N; ++
n) {
3886 Int16 const * p = mLeafNodes[
n]->buffer().data();
3890 count +=
Index32(sEdgeGroupTable[(SIGNS & *p)][0]);
3899 Int16LeafNodeType *
const *
const mLeafNodes;
3904 template<
typename Po
intIndexLeafNode>
3905 struct AdaptiveLeafNodePointCount
3907 using Int16LeafNodeType = tree::LeafNode<Int16, PointIndexLeafNode::LOG2DIM>;
3909 AdaptiveLeafNodePointCount(
const std::vector<PointIndexLeafNode*>& pointIndexNodes,
3910 const std::vector<Int16LeafNodeType*>& signDataNodes,
3911 std::unique_ptr<
Index32[]>& leafNodeCount)
3912 : mPointIndexNodes(pointIndexNodes.data())
3913 , mSignDataNodes(signDataNodes.data())
3914 , mData(leafNodeCount.
get()) {}
3916 void operator()(
const tbb::blocked_range<size_t>& range)
const
3920 for (
size_t n = range.begin(), N = range.end(); n !=
N; ++
n) {
3922 const PointIndexLeafNode& node = *mPointIndexNodes[
n];
3923 const Int16LeafNodeType& signNode = *mSignDataNodes[
n];
3927 std::set<IndexType> uniqueRegions;
3929 for (
typename PointIndexLeafNode::ValueOnCIter it = node.cbeginValueOn(); it; ++it) {
3934 count += size_t(sEdgeGroupTable[(SIGNS & signNode.getValue(it.pos()))][0]);
3936 uniqueRegions.insert(
id);
3940 mData[
n] =
Index32(count + uniqueRegions.size());
3945 PointIndexLeafNode
const *
const *
const mPointIndexNodes;
3946 Int16LeafNodeType
const *
const *
const mSignDataNodes;
3951 template<
typename Po
intIndexLeafNode>
3954 using Int16LeafNodeType = tree::LeafNode<Int16, PointIndexLeafNode::LOG2DIM>;
3956 MapPoints(
const std::vector<PointIndexLeafNode*>& pointIndexNodes,
3957 const std::vector<Int16LeafNodeType*>& signDataNodes,
3958 std::unique_ptr<
Index32[]>& leafNodeCount)
3959 : mPointIndexNodes(pointIndexNodes.data())
3960 , mSignDataNodes(signDataNodes.data())
3961 , mData(leafNodeCount.
get()) {}
3963 void operator()(
const tbb::blocked_range<size_t>& range)
const {
3965 for (
size_t n = range.begin(), N = range.end(); n !=
N; ++
n) {
3967 const Int16LeafNodeType& signNode = *mSignDataNodes[
n];
3968 PointIndexLeafNode& indexNode = *mPointIndexNodes[
n];
3972 for (
auto it = indexNode.beginValueOn(); it; ++it) {
3973 const Index pos = it.pos();
3974 indexNode.setValueOnly(pos, pointOffset);
3975 const int signs = SIGNS &
int(signNode.getValue(pos));
3976 pointOffset +=
Index32(sEdgeGroupTable[signs][0]);
3982 PointIndexLeafNode *
const *
const mPointIndexNodes;
3983 Int16LeafNodeType
const *
const *
const mSignDataNodes;
3988 template<
typename TreeType,
typename PrimBuilder>
3989 struct ComputePolygons
3992 using Int16LeafNodeType =
typename Int16TreeType::LeafNodeType;
3995 using Index32LeafNodeType =
typename Index32TreeType::LeafNodeType;
3998 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
3999 const Int16TreeType& signFlagsTree,
4000 const Index32TreeType& idxTree,
4002 bool invertSurfaceOrientation);
4004 void setRefSignTree(
const Int16TreeType *
r) { mRefSignFlagsTree =
r; }
4006 void operator()(
const tbb::blocked_range<size_t>&)
const;
4009 Int16LeafNodeType *
const *
const mSignFlagsLeafNodes;
4010 Int16TreeType
const *
const mSignFlagsTree;
4011 Int16TreeType
const * mRefSignFlagsTree;
4012 Index32TreeType
const *
const mIndexTree;
4014 bool const mInvertSurfaceOrientation;
4018 template<
typename TreeType,
typename PrimBuilder>
4019 ComputePolygons<TreeType, PrimBuilder>::ComputePolygons(
4020 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
4021 const Int16TreeType& signFlagsTree,
4022 const Index32TreeType& idxTree,
4024 bool invertSurfaceOrientation)
4025 : mSignFlagsLeafNodes(signFlagsLeafNodes.data())
4026 , mSignFlagsTree(&signFlagsTree)
4027 , mRefSignFlagsTree(nullptr)
4028 , mIndexTree(&idxTree)
4029 , mPolygonPoolList(&polygons)
4030 , mInvertSurfaceOrientation(invertSurfaceOrientation)
4034 template<
typename InputTreeType,
typename PrimBuilder>
4036 ComputePolygons<InputTreeType, PrimBuilder>::operator()(
const tbb::blocked_range<size_t>& range)
const
4038 using Int16ValueAccessor = tree::ValueAccessor<const Int16TreeType>;
4039 Int16ValueAccessor signAcc(*mSignFlagsTree);
4041 tree::ValueAccessor<const Index32TreeType> idxAcc(*mIndexTree);
4043 const bool invertSurfaceOrientation = mInvertSurfaceOrientation;
4050 std::unique_ptr<Int16ValueAccessor> refSignAcc;
4051 if (mRefSignFlagsTree) refSignAcc.reset(
new Int16ValueAccessor(*mRefSignFlagsTree));
4053 for (
size_t n = range.begin(); n != range.end(); ++
n) {
4055 const Int16LeafNodeType& node = *mSignFlagsLeafNodes[
n];
4056 origin = node.origin();
4060 typename Int16LeafNodeType::ValueOnCIter iter = node.cbeginValueOn();
4061 for (; iter; ++iter) {
4062 if (iter.getValue() & XEDGE) ++edgeCount;
4063 if (iter.getValue() & YEDGE) ++edgeCount;
4064 if (iter.getValue() & ZEDGE) ++edgeCount;
4067 if (edgeCount == 0)
continue;
4069 mesher.init(edgeCount, (*mPolygonPoolList)[n]);
4071 const Int16LeafNodeType *signleafPt = signAcc.probeConstLeaf(origin);
4072 const Index32LeafNodeType *idxLeafPt = idxAcc.probeConstLeaf(origin);
4074 if (!signleafPt || !idxLeafPt)
continue;
4076 const Int16LeafNodeType *refSignLeafPt =
nullptr;
4077 if (refSignAcc) refSignLeafPt = refSignAcc->probeConstLeaf(origin);
4081 for (iter = node.cbeginValueOn(); iter; ++iter) {
4082 ijk = iter.getCoord();
4084 const Int16 flags = iter.getValue();
4085 if (!(flags & 0xE00))
continue;
4088 if (refSignLeafPt) {
4089 refFlags = refSignLeafPt->getValue(iter.pos());
4092 const uint8_t cell = uint8_t(SIGNS & flags);
4094 if (sEdgeGroupTable[cell][0] > 1) {
4095 offsets[0] = (sEdgeGroupTable[cell][1] - 1);
4096 offsets[1] = (sEdgeGroupTable[cell][9] - 1);
4097 offsets[2] = (sEdgeGroupTable[cell][4] - 1);
4103 if (ijk[0] > origin[0] && ijk[1] > origin[1] && ijk[2] > origin[2]) {
4104 constructPolygons(invertSurfaceOrientation,
4105 flags, refFlags, offsets, ijk, *signleafPt, *idxLeafPt, mesher);
4107 constructPolygons(invertSurfaceOrientation,
4108 flags, refFlags, offsets, ijk, signAcc, idxAcc, mesher);
4121 template<
typename T>
4124 CopyArray(T * outputArray,
const T * inputArray,
size_t outputOffset = 0)
4125 : mOutputArray(outputArray), mInputArray(inputArray), mOutputOffset(outputOffset)
4129 void operator()(
const tbb::blocked_range<size_t>& inputArrayRange)
const
4131 const size_t offset = mOutputOffset;
4132 for (
size_t n = inputArrayRange.begin(), N = inputArrayRange.end(); n <
N; ++
n) {
4133 mOutputArray[offset +
n] = mInputArray[
n];
4138 T *
const mOutputArray;
4139 T const *
const mInputArray;
4140 size_t const mOutputOffset;
4144 struct FlagAndCountQuadsToSubdivide
4147 const std::vector<uint8_t>& pointFlags,
4149 std::unique_ptr<
unsigned[]>& numQuadsToDivide)
4150 : mPolygonPoolList(&polygons)
4151 , mPointFlags(pointFlags.data())
4152 , mPoints(points.
get())
4153 , mNumQuadsToDivide(numQuadsToDivide.
get())
4157 void operator()(
const tbb::blocked_range<size_t>& range)
const
4159 for (
size_t n = range.begin(), N = range.end(); n <
N; ++
n) {
4161 PolygonPool& polygons = (*mPolygonPoolList)[
n];
4166 for (
size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4168 char& flags = polygons.quadFlags(i);
4172 Vec4I& quad = polygons.quad(i);
4174 const bool edgePoly = mPointFlags[quad[0]] || mPointFlags[quad[1]]
4175 || mPointFlags[quad[2]] || mPointFlags[quad[3]];
4177 if (!edgePoly)
continue;
4179 const Vec3s& p0 = mPoints[quad[0]];
4180 const Vec3s& p1 = mPoints[quad[1]];
4181 const Vec3s& p2 = mPoints[quad[2]];
4182 const Vec3s& p3 = mPoints[quad[3]];
4184 if (!isPlanarQuad(p0, p1, p2, p3, 1e-6
f)) {
4191 mNumQuadsToDivide[
n] = count;
4197 uint8_t
const *
const mPointFlags;
4198 Vec3s const *
const mPoints;
4199 unsigned *
const mNumQuadsToDivide;
4203 struct SubdivideQuads
4209 std::unique_ptr<
unsigned[]>& numQuadsToDivide,
4210 std::unique_ptr<
unsigned[]>& centroidOffsets)
4211 : mPolygonPoolList(&polygons)
4212 , mPoints(points.
get())
4213 , mCentroids(centroids.
get())
4214 , mNumQuadsToDivide(numQuadsToDivide.
get())
4215 , mCentroidOffsets(centroidOffsets.
get())
4220 void operator()(
const tbb::blocked_range<size_t>& range)
const
4222 for (
size_t n = range.begin(), N = range.end(); n <
N; ++
n) {
4224 PolygonPool& polygons = (*mPolygonPoolList)[
n];
4226 const size_t nonplanarCount = size_t(mNumQuadsToDivide[n]);
4228 if (nonplanarCount > 0) {
4230 PolygonPool tmpPolygons;
4231 tmpPolygons.resetQuads(polygons.numQuads() - nonplanarCount);
4232 tmpPolygons.resetTriangles(polygons.numTriangles() + size_t(4) * nonplanarCount);
4234 size_t offset = mCentroidOffsets[
n];
4236 size_t triangleIdx = 0;
4238 for (
size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4240 const char quadFlags = polygons.quadFlags(i);
4243 unsigned newPointIdx = unsigned(offset +
mPointCount);
4247 mCentroids[
offset] = (mPoints[quad[0]] + mPoints[quad[1]] +
4248 mPoints[quad[2]] + mPoints[quad[3]]) * 0.25
f;
4253 Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
4255 triangle[0] = quad[0];
4256 triangle[1] = newPointIdx;
4257 triangle[2] = quad[3];
4259 tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
4265 Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
4267 triangle[0] = quad[0];
4268 triangle[1] = quad[1];
4269 triangle[2] = newPointIdx;
4271 tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
4277 Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
4279 triangle[0] = quad[1];
4280 triangle[1] = quad[2];
4281 triangle[2] = newPointIdx;
4283 tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
4290 Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
4292 triangle[0] = quad[2];
4293 triangle[1] = quad[3];
4294 triangle[2] = newPointIdx;
4296 tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
4305 for (
size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
4306 tmpPolygons.triangle(triangleIdx) = polygons.triangle(i);
4307 tmpPolygons.triangleFlags(triangleIdx) = polygons.triangleFlags(i);
4312 for (
size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4316 tmpPolygons.quad(quadIdx) = quad;
4317 tmpPolygons.quadFlags(quadIdx) = polygons.quadFlags(i);
4322 polygons.copy(tmpPolygons);
4329 Vec3s const *
const mPoints;
4330 Vec3s *
const mCentroids;
4331 unsigned *
const mNumQuadsToDivide;
4332 unsigned *
const mCentroidOffsets;
4338 subdivideNonplanarSeamLineQuads(
4340 size_t polygonPoolListSize,
4342 size_t& pointListSize,
4343 std::vector<uint8_t>& pointFlags)
4345 const tbb::blocked_range<size_t> polygonPoolListRange(0, polygonPoolListSize);
4347 std::unique_ptr<unsigned[]> numQuadsToDivide(
new unsigned[polygonPoolListSize]);
4350 FlagAndCountQuadsToSubdivide(polygonPoolList, pointFlags, pointList, numQuadsToDivide));
4352 std::unique_ptr<unsigned[]> centroidOffsets(
new unsigned[polygonPoolListSize]);
4354 size_t centroidCount = 0;
4358 for (
size_t n = 0, N = polygonPoolListSize; n <
N; ++
n) {
4359 centroidOffsets[
n] = sum;
4360 sum += numQuadsToDivide[
n];
4362 centroidCount = size_t(sum);
4365 std::unique_ptr<Vec3s[]> centroidList(
new Vec3s[centroidCount]);
4368 SubdivideQuads(polygonPoolList, pointList, pointListSize,
4369 centroidList, numQuadsToDivide, centroidOffsets));
4371 if (centroidCount > 0) {
4373 const size_t newPointListSize = centroidCount + pointListSize;
4375 std::unique_ptr<openvdb::Vec3s[]> newPointList(
new openvdb::Vec3s[newPointListSize]);
4378 CopyArray<Vec3s>(newPointList.get(), pointList.get()));
4380 tbb::parallel_for(tbb::blocked_range<size_t>(0, newPointListSize - pointListSize),
4381 CopyArray<Vec3s>(newPointList.get(), centroidList.get(), pointListSize));
4383 pointListSize = newPointListSize;
4384 pointList.swap(newPointList);
4385 pointFlags.resize(pointListSize, 0);
4390 struct ReviseSeamLineFlags
4393 const std::vector<uint8_t>& pointFlags)
4394 : mPolygonPoolList(&polygons)
4395 , mPointFlags(pointFlags.data())
4399 void operator()(
const tbb::blocked_range<size_t>& range)
const
4401 for (
size_t n = range.begin(), N = range.end(); n <
N; ++
n) {
4403 PolygonPool& polygons = (*mPolygonPoolList)[
n];
4405 for (
size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4407 char& flags = polygons.quadFlags(i);
4409 if (flags & POLYFLAG_FRACTURE_SEAM) {
4413 const bool hasSeamLinePoint =
4414 mPointFlags[verts[0]] || mPointFlags[verts[1]] ||
4415 mPointFlags[verts[2]] || mPointFlags[verts[3]];
4417 if (!hasSeamLinePoint) {
4418 flags &= ~POLYFLAG_FRACTURE_SEAM;
4423 for (
size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
4425 char& flags = polygons.triangleFlags(i);
4427 if (flags & POLYFLAG_FRACTURE_SEAM) {
4431 const bool hasSeamLinePoint =
4432 mPointFlags[verts[0]] || mPointFlags[verts[1]] || mPointFlags[verts[2]];
4434 if (!hasSeamLinePoint) {
4435 flags &= ~POLYFLAG_FRACTURE_SEAM;
4446 uint8_t
const *
const mPointFlags;
4451 reviseSeamLineFlags(
PolygonPoolList& polygonPoolList,
size_t polygonPoolListSize,
4452 std::vector<uint8_t>& pointFlags)
4455 ReviseSeamLineFlags(polygonPoolList, pointFlags));
4462 template<
typename InputTreeType>
4463 struct MaskDisorientedTrianglePoints
4465 MaskDisorientedTrianglePoints(
const InputTreeType& inputTree,
const PolygonPoolList& polygons,
4466 const PointList& pointList, std::unique_ptr<uint8_t[]>& pointMask,
4467 const math::Transform& transform,
bool invertSurfaceOrientation)
4468 : mInputTree(&inputTree)
4469 , mPolygonPoolList(&polygons)
4470 , mPointList(&pointList)
4471 , mPointMask(pointMask.
get())
4472 , mTransform(transform)
4473 , mInvertSurfaceOrientation(invertSurfaceOrientation)
4477 void operator()(
const tbb::blocked_range<size_t>& range)
const
4481 tree::ValueAccessor<const InputTreeType> inputAcc(*mInputTree);
4482 Vec3s centroid, normal;
4485 const bool invertGradientDir = mInvertSurfaceOrientation || isBoolValue<ValueType>();
4487 for (
size_t n = range.begin(), N = range.end(); n <
N; ++
n) {
4489 const PolygonPool& polygons = (*mPolygonPoolList)[
n];
4491 for (
size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
4493 const Vec3I& verts = polygons.triangle(i);
4495 const Vec3s& v0 = (*mPointList)[verts[0]];
4496 const Vec3s& v1 = (*mPointList)[verts[1]];
4497 const Vec3s& v2 = (*mPointList)[verts[2]];
4499 normal = (v2 -
v0).
cross((v1 - v0));
4502 centroid = (v0 + v1 +
v2) * (1.0
f / 3.0
f);
4503 ijk = mTransform.worldToIndexCellCentered(centroid);
4508 if (invertGradientDir) {
4513 if (dir.dot(normal) < -0.5f) {
4518 mPointMask[verts[0]] = 1;
4519 mPointMask[verts[1]] = 1;
4520 mPointMask[verts[2]] = 1;
4529 InputTreeType
const *
const mInputTree;
4532 uint8_t *
const mPointMask;
4533 math::Transform
const mTransform;
4534 bool const mInvertSurfaceOrientation;
4538 template<
typename InputTree>
4540 relaxDisorientedTriangles(
4541 bool invertSurfaceOrientation,
4542 const InputTree& inputTree,
4543 const math::Transform& transform,
4545 size_t polygonPoolListSize,
4547 const size_t pointListSize)
4549 const tbb::blocked_range<size_t> polygonPoolListRange(0, polygonPoolListSize);
4551 std::unique_ptr<uint8_t[]> pointMask(
new uint8_t[pointListSize]);
4552 fillArray(pointMask.get(), uint8_t(0), pointListSize);
4555 MaskDisorientedTrianglePoints<InputTree>(
4556 inputTree, polygonPoolList, pointList, pointMask, transform, invertSurfaceOrientation));
4558 std::unique_ptr<uint8_t[]> pointUpdates(
new uint8_t[pointListSize]);
4559 fillArray(pointUpdates.get(), uint8_t(0), pointListSize);
4561 std::unique_ptr<Vec3s[]> newPoints(
new Vec3s[pointListSize]);
4562 fillArray(newPoints.get(),
Vec3s(0.0
f, 0.0
f, 0.0
f), pointListSize);
4564 for (
size_t n = 0, N = polygonPoolListSize; n <
N; ++
n) {
4566 PolygonPool& polygons = polygonPoolList[
n];
4568 for (
size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4571 for (
int v = 0; v < 4; ++
v) {
4573 const unsigned pointIdx = verts[
v];
4575 if (pointMask[pointIdx] == 1) {
4577 newPoints[pointIdx] +=
4578 pointList[verts[0]] + pointList[verts[1]] +
4579 pointList[verts[2]] + pointList[verts[3]];
4581 pointUpdates[pointIdx] = uint8_t(pointUpdates[pointIdx] + 4);
4586 for (
size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
4589 for (
int v = 0; v < 3; ++
v) {
4591 const unsigned pointIdx = verts[
v];
4593 if (pointMask[pointIdx] == 1) {
4594 newPoints[pointIdx] +=
4595 pointList[verts[0]] + pointList[verts[1]] + pointList[verts[2]];
4597 pointUpdates[pointIdx] = uint8_t(pointUpdates[pointIdx] + 3);
4603 for (
size_t n = 0, N = pointListSize; n <
N; ++
n) {
4604 if (pointUpdates[n] > 0) {
4605 const double weight = 1.0 / double(pointUpdates[n]);
4606 pointList[
n] = newPoints[
n] *
float(weight);
4612 template<
typename Gr
idType>
4616 std::vector<Vec3s>& points,
4617 std::vector<Vec3I>& triangles,
4618 std::vector<Vec4I>&
quads,
4621 bool relaxDisorientedTriangles)
4624 "volume to mesh conversion is supported only for scalar grids");
4626 VolumeToMesh mesher(isovalue, adaptivity, relaxDisorientedTriangles);
4631 points.resize(mesher.pointListSize());
4634 volume_to_mesh_internal::PointListCopy ptnCpy(mesher.pointList(),
points);
4636 mesher.pointList().reset(
nullptr);
4642 size_t numQuads = 0, numTriangles = 0;
4643 for (
size_t n = 0, N = mesher.polygonPoolListSize(); n <
N; ++
n) {
4644 openvdb::tools::PolygonPool& polygons = polygonPoolList[
n];
4645 numTriangles += polygons.numTriangles();
4646 numQuads += polygons.numQuads();
4650 triangles.resize(numTriangles);
4652 quads.resize(numQuads);
4656 size_t qIdx = 0, tIdx = 0;
4657 for (
size_t n = 0, N = mesher.polygonPoolListSize(); n <
N; ++
n) {
4658 openvdb::tools::PolygonPool& polygons = polygonPoolList[
n];
4660 for (
size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4661 quads[qIdx++] = polygons.quad(i);
4664 for (
size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
4665 triangles[tIdx++] = polygons.triangle(i);
4679 PolygonPool::PolygonPool()
4683 , mTriangles(nullptr)
4684 , mQuadFlags(nullptr)
4685 , mTriangleFlags(nullptr)
4692 : mNumQuads(numQuads)
4693 , mNumTriangles(numTriangles)
4694 , mQuads(new openvdb::
Vec4I[mNumQuads])
4695 , mTriangles(new openvdb::
Vec3I[mNumTriangles])
4696 , mQuadFlags(new char[mNumQuads])
4697 , mTriangleFlags(new char[mNumTriangles])
4708 for (
size_t i = 0; i < mNumQuads; ++i) {
4709 mQuads[i] = rhs.mQuads[i];
4710 mQuadFlags[i] = rhs.mQuadFlags[i];
4713 for (
size_t i = 0; i < mNumTriangles; ++i) {
4714 mTriangles[i] = rhs.mTriangles[i];
4715 mTriangleFlags[i] = rhs.mTriangleFlags[i];
4725 mQuadFlags.reset(
new char[mNumQuads]);
4733 mQuads.reset(
nullptr);
4734 mQuadFlags.reset(
nullptr);
4741 mNumTriangles =
size;
4743 mTriangleFlags.reset(
new char[mNumTriangles]);
4751 mTriangles.reset(
nullptr);
4752 mTriangleFlags.reset(
nullptr);
4759 if (!(n < mNumQuads))
return false;
4764 mQuads.reset(
nullptr);
4768 std::unique_ptr<char[]>
flags(
new char[n]);
4770 for (
size_t i = 0; i <
n; ++i) {
4771 quads[i] = mQuads[i];
4772 flags[i] = mQuadFlags[i];
4776 mQuadFlags.swap(flags);
4788 if (!(n < mNumTriangles))
return false;
4793 mTriangles.reset(
nullptr);
4796 std::unique_ptr<openvdb::Vec3I[]> triangles(
new openvdb::Vec3I[n]);
4797 std::unique_ptr<char[]>
flags(
new char[n]);
4799 for (
size_t i = 0; i <
n; ++i) {
4800 triangles[i] = mTriangles[i];
4801 flags[i] = mTriangleFlags[i];
4804 mTriangles.swap(triangles);
4805 mTriangleFlags.swap(flags);
4822 , mSeamPointListSize(0)
4823 , mPolygonPoolListSize(0)
4824 , mIsovalue(isovalue)
4825 , mPrimAdaptivity(adaptivity)
4826 , mSecAdaptivity(0.0)
4828 , mSurfaceMaskGrid(
GridBase::ConstPtr())
4829 , mAdaptivityGrid(
GridBase::ConstPtr())
4830 , mAdaptivityMaskTree(
TreeBase::ConstPtr())
4833 , mInvertSurfaceMask(false)
4834 , mRelaxDisorientedTriangles(relaxDisorientedTriangles)
4835 , mQuantizedSeamPoints(nullptr)
4845 mSecAdaptivity = secAdaptivity;
4850 mSeamPointListSize = 0;
4851 mQuantizedSeamPoints.reset(
nullptr);
4858 mSurfaceMaskGrid =
mask;
4859 mInvertSurfaceMask = invertMask;
4866 mAdaptivityGrid = grid;
4873 mAdaptivityMaskTree = tree;
4877 template<
typename InputGr
idType>
4883 using InputTreeType =
typename InputGridType::TreeType;
4884 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
4891 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
4893 using Int16LeafNodeType =
typename Int16TreeType::LeafNodeType;
4895 using Index32LeafNodeType =
typename Index32TreeType::LeafNodeType;
4900 mPolygonPoolListSize = 0;
4902 mPointFlags.clear();
4907 const InputValueType isovalue = InputValueType(mIsovalue);
4908 const float adaptivityThreshold =
float(mPrimAdaptivity);
4909 const bool adaptive = mPrimAdaptivity > 1e-7 || mSecAdaptivity > 1e-7;
4915 const bool invertSurfaceOrientation = (!volume_to_mesh_internal::isBoolValue<InputValueType>()
4920 const InputTreeType& inputTree = inputGrid.tree();
4922 BoolTreeType intersectionTree(
false), adaptivityMask(
false);
4924 if (mAdaptivityMaskTree && mAdaptivityMaskTree->type() == BoolTreeType::treeType()) {
4925 const BoolTreeType *refAdaptivityMask=
4926 static_cast<const BoolTreeType*
>(mAdaptivityMaskTree.get());
4927 adaptivityMask.topologyUnion(*refAdaptivityMask);
4930 Int16TreeType signFlagsTree(0);
4936 volume_to_mesh_internal::identifySurfaceIntersectingVoxels(
4937 intersectionTree, inputTree, isovalue);
4939 volume_to_mesh_internal::applySurfaceMask(intersectionTree, adaptivityMask,
4940 inputGrid, mSurfaceMaskGrid, mInvertSurfaceMask, isovalue);
4942 if (intersectionTree.empty())
return;
4944 volume_to_mesh_internal::computeAuxiliaryData(
4945 signFlagsTree, pointIndexTree, intersectionTree, inputTree, isovalue);
4947 intersectionTree.clear();
4949 std::vector<Index32LeafNodeType*> pointIndexLeafNodes;
4950 pointIndexTree.getNodes(pointIndexLeafNodes);
4952 std::vector<Int16LeafNodeType*> signFlagsLeafNodes;
4953 signFlagsTree.getNodes(signFlagsLeafNodes);
4955 const tbb::blocked_range<size_t> auxiliaryLeafNodeRange(0, signFlagsLeafNodes.size());
4960 Int16TreeType* refSignFlagsTree =
nullptr;
4961 Index32TreeType* refPointIndexTree =
nullptr;
4962 InputTreeType
const* refInputTree =
nullptr;
4964 if (mRefGrid && mRefGrid->type() == InputGridType::gridType()) {
4966 const InputGridType* refGrid =
static_cast<const InputGridType*
>(mRefGrid.get());
4967 refInputTree = &refGrid->tree();
4969 if (!mRefSignTree && !mRefIdxTree) {
4973 typename Int16TreeType::Ptr refSignFlagsTreePt(
new Int16TreeType(0));
4974 typename Index32TreeType::Ptr refPointIndexTreePt(
4977 BoolTreeType refIntersectionTree(
false);
4979 volume_to_mesh_internal::identifySurfaceIntersectingVoxels(
4980 refIntersectionTree, *refInputTree, isovalue);
4982 volume_to_mesh_internal::computeAuxiliaryData(*refSignFlagsTreePt,
4983 *refPointIndexTreePt, refIntersectionTree, *refInputTree, isovalue);
4985 mRefSignTree = refSignFlagsTreePt;
4986 mRefIdxTree = refPointIndexTreePt;
4989 if (mRefSignTree && mRefIdxTree) {
4993 refSignFlagsTree =
static_cast<Int16TreeType*
>(mRefSignTree.get());
4994 refPointIndexTree =
static_cast<Index32TreeType*
>(mRefIdxTree.get());
4998 if (refSignFlagsTree && refPointIndexTree) {
5002 volume_to_mesh_internal::markSeamLineData(signFlagsTree, *refSignFlagsTree);
5004 if (mSeamPointListSize == 0) {
5008 std::vector<Int16LeafNodeType*> refSignFlagsLeafNodes;
5009 refSignFlagsTree->getNodes(refSignFlagsLeafNodes);
5011 std::unique_ptr<Index32[]> leafNodeOffsets(
5012 new Index32[refSignFlagsLeafNodes.size()]);
5015 volume_to_mesh_internal::LeafNodePointCount<Int16LeafNodeType::LOG2DIM>(
5016 refSignFlagsLeafNodes, leafNodeOffsets));
5020 for (
size_t n = 0, N = refSignFlagsLeafNodes.size(); n <
N; ++
n) {
5021 const Index32 tmp = leafNodeOffsets[
n];
5022 leafNodeOffsets[
n] = count;
5025 mSeamPointListSize = size_t(count);
5028 if (mSeamPointListSize != 0) {
5030 mQuantizedSeamPoints.reset(
new uint32_t[mSeamPointListSize]);
5032 std::memset(mQuantizedSeamPoints.get(), 0,
sizeof(uint32_t) * mSeamPointListSize);
5034 std::vector<Index32LeafNodeType*> refPointIndexLeafNodes;
5035 refPointIndexTree->getNodes(refPointIndexLeafNodes);
5038 volume_to_mesh_internal::MapPoints<Index32LeafNodeType>(
5039 refPointIndexLeafNodes, refSignFlagsLeafNodes, leafNodeOffsets));
5043 if (mSeamPointListSize != 0) {
5046 volume_to_mesh_internal::SeamLineWeights<InputTreeType>(
5047 signFlagsLeafNodes, inputTree, *refPointIndexTree, *refSignFlagsTree,
5048 mQuantizedSeamPoints.get(), isovalue));
5053 const bool referenceMeshing = refSignFlagsTree && refPointIndexTree && refInputTree;
5058 std::unique_ptr<Index32[]> leafNodeOffsets(
new Index32[signFlagsLeafNodes.size()]);
5061 volume_to_mesh_internal::MergeVoxelRegions<InputGridType> mergeOp(
5062 inputGrid, pointIndexTree, pointIndexLeafNodes, signFlagsLeafNodes,
5063 isovalue, adaptivityThreshold, invertSurfaceOrientation);
5065 if (mAdaptivityGrid && mAdaptivityGrid->type() == FloatGridType::gridType()) {
5066 const FloatGridType* adaptivityGrid =
5067 static_cast<const FloatGridType*
>(mAdaptivityGrid.get());
5068 mergeOp.setSpatialAdaptivity(*adaptivityGrid);
5071 if (!adaptivityMask.empty()) {
5072 mergeOp.setAdaptivityMask(adaptivityMask);
5075 if (referenceMeshing) {
5076 mergeOp.setRefSignFlagsData(*refSignFlagsTree,
float(mSecAdaptivity));
5081 volume_to_mesh_internal::AdaptiveLeafNodePointCount<Index32LeafNodeType>
5082 op(pointIndexLeafNodes, signFlagsLeafNodes, leafNodeOffsets);
5088 volume_to_mesh_internal::LeafNodePointCount<Int16LeafNodeType::LOG2DIM>
5089 op(signFlagsLeafNodes, leafNodeOffsets);
5097 for (
size_t n = 0, N = signFlagsLeafNodes.size(); n <
N; ++
n) {
5098 const Index32 tmp = leafNodeOffsets[
n];
5103 mPointListSize = size_t(pointCount);
5105 mPointFlags.clear();
5112 volume_to_mesh_internal::ComputePoints<InputTreeType>
5113 op(mPoints.get(), inputTree, pointIndexLeafNodes,
5114 signFlagsLeafNodes, leafNodeOffsets,
transform, mIsovalue);
5116 if (referenceMeshing) {
5117 mPointFlags.resize(mPointListSize);
5118 op.setRefData(*refInputTree, *refPointIndexTree, *refSignFlagsTree,
5119 mQuantizedSeamPoints.get(), mPointFlags.data());
5128 mPolygonPoolListSize = signFlagsLeafNodes.size();
5129 mPolygons.reset(
new PolygonPool[mPolygonPoolListSize]);
5133 using PrimBuilder = volume_to_mesh_internal::AdaptivePrimBuilder;
5135 volume_to_mesh_internal::ComputePolygons<Int16TreeType, PrimBuilder>
5136 op(signFlagsLeafNodes, signFlagsTree, pointIndexTree,
5137 mPolygons, invertSurfaceOrientation);
5139 if (referenceMeshing) {
5140 op.setRefSignTree(refSignFlagsTree);
5147 using PrimBuilder = volume_to_mesh_internal::UniformPrimBuilder;
5149 volume_to_mesh_internal::ComputePolygons<Int16TreeType, PrimBuilder>
5150 op(signFlagsLeafNodes, signFlagsTree, pointIndexTree,
5151 mPolygons, invertSurfaceOrientation);
5153 if (referenceMeshing) {
5154 op.setRefSignTree(refSignFlagsTree);
5161 signFlagsTree.clear();
5162 pointIndexTree.clear();
5165 if (adaptive && mRelaxDisorientedTriangles) {
5166 volume_to_mesh_internal::relaxDisorientedTriangles(invertSurfaceOrientation,
5167 inputTree, transform, mPolygons, mPolygonPoolListSize, mPoints, mPointListSize);
5171 if (referenceMeshing) {
5172 volume_to_mesh_internal::subdivideNonplanarSeamLineQuads(
5173 mPolygons, mPolygonPoolListSize, mPoints, mPointListSize, mPointFlags);
5175 volume_to_mesh_internal::reviseSeamLineFlags(mPolygons, mPolygonPoolListSize, mPointFlags);
5184 template<
typename Gr
idType>
5187 std::vector<Vec3s>& points,
5188 std::vector<Vec3I>& triangles,
5189 std::vector<Vec4I>& quads,
5192 bool relaxDisorientedTriangles)
5194 volume_to_mesh_internal::doVolumeToMesh(grid, points, triangles, quads,
5195 isovalue, adaptivity, relaxDisorientedTriangles);
5198 template<
typename Gr
idType>
5201 std::vector<Vec3s>& points,
5202 std::vector<Vec4I>& quads,
5205 std::vector<Vec3I> triangles;
5206 volumeToMesh(grid, points, triangles, quads, isovalue, 0.0,
true);
5215 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
5217 #ifdef OPENVDB_INSTANTIATE_VOLUMETOMESH
5221 #define _FUNCTION(TreeT) \
5222 void volumeToMesh(const Grid<TreeT>&, std::vector<Vec3s>&, std::vector<Vec4I>&, double)
5226 #define _FUNCTION(TreeT) \
5227 void volumeToMesh(const Grid<TreeT>&, std::vector<Vec3s>&, std::vector<Vec3I>&, std::vector<Vec4I>&, double, double, bool)
5231 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
5238 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
void parallel_for(int64_t start, int64_t end, std::function< void(int64_t index)> &&task, parallel_options opt=parallel_options(0, Split_Y, 1))
GU_API exint quads(GU_Detail *gdp, UT_Array< GA_OffsetArray > &rings, UT_Array< GA_OffsetArray > &originalRings, GA_PrimitiveGroup *patchgroup, GA_PrimitiveGroup *loopgroup, bool smooth, fpreal smoothstrength, bool edgeloop, fpreal edgeloopPercentage)
typedef int(APIENTRYP RE_PFNGLXSWAPINTERVALSGIPROC)(int)
__hostdev__ uint64_t pointCount() const
SharedPtr< TreeBase > Ptr
GLdouble GLdouble GLint GLint const GLdouble * points
constexpr Index32 INVALID_IDX
GridType
List of types that are currently supported by NanoVDB.
GLsizei const GLfloat * value
**And then you can **find out if it s done
GLdouble GLdouble GLdouble z
void reverse(I begin, I end)
GLuint GLsizei GLsizei * length
#define OPENVDB_USE_VERSION_NAMESPACE
ImageBuf OIIO_API min(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
**But if you need a result
GLfloat GLfloat GLfloat v2
GLuint GLsizei const GLuint const GLintptr * offsets
SIM_DerVector3 normalize() const
Abstract base class for typed grids.
Mat3 transpose() const
returns transpose of this
SYS_FORCE_INLINE const_iterator end() const
std::vector< Index > IndexArray
IMATH_NAMESPACE::V2f float
math::Vec4< Index32 > Vec4I
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) ...
ValueAccessors are designed to help accelerate accesses into the OpenVDB Tree structures by storing c...
auto get(const UT_ARTIterator< T > &it) -> decltype(it.key())
GA_API const UT_StringHolder transform
Container class that associates a tree with a transform and metadata.
#define OPENVDB_NUMERIC_TREE_INSTANTIATE(Function)
Base class for typed trees.
__hostdev__ T dot(const Vec3T &v) const
GLenum GLsizei GLsizei GLint * values
SharedPtr< const GridBase > ConstPtr
math::Vec3< Index32 > Vec3I
ImageBuf OIIO_API max(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
GA_API const UT_StringHolder N
SharedPtr< const TreeBase > ConstPtr
GLubyte GLubyte GLubyte GLubyte w
IMATH_INTERNAL_NAMESPACE_HEADER_ENTER IMATH_HOSTDEVICE constexpr T abs(T a) IMATH_NOEXCEPT
void OIIO_UTIL_API split(string_view str, std::vector< string_view > &result, string_view sep=string_view(), int maxsplit=-1)
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
SIM_DerVector3 cross(const SIM_DerVector3 &lhs, const SIM_DerVector3 &rhs)
#define OPENVDB_THROW(exception, message)