10 #ifndef OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
11 #define OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
18 #include <tbb/blocked_range.h>
19 #include <tbb/parallel_for.h>
20 #include <tbb/parallel_reduce.h>
21 #include <tbb/task_arena.h>
27 #include <type_traits>
51 template<
typename Gr
idType>
55 std::vector<Vec3s>&
points,
56 std::vector<Vec4I>&
quads,
57 double isovalue = 0.0);
72 template<
typename Gr
idType>
76 std::vector<Vec3s>&
points,
77 std::vector<Vec3I>& triangles,
78 std::vector<Vec4I>&
quads,
79 double isovalue = 0.0,
80 double adaptivity = 0.0,
81 bool relaxDisorientedTriangles =
true);
110 const size_t&
numQuads()
const {
return mNumQuads; }
134 inline bool trimQuads(
const size_t n,
bool reallocate =
false);
135 inline bool trimTrinagles(
const size_t n,
bool reallocate =
false);
141 size_t mNumQuads, mNumTriangles;
142 std::unique_ptr<openvdb::Vec4I[]> mQuads;
143 std::unique_ptr<openvdb::Vec3I[]> mTriangles;
144 std::unique_ptr<char[]> mQuadFlags, mTriangleFlags;
166 VolumeToMesh(
double isovalue = 0,
double adaptivity = 0,
bool relaxDisorientedTriangles =
true);
182 const std::vector<uint8_t>&
pointFlags()
const {
return mPointFlags; }
191 template<
typename InputGr
idType>
245 size_t mPointListSize, mSeamPointListSize, mPolygonPoolListSize;
246 double mIsovalue, mPrimAdaptivity, mSecAdaptivity;
253 bool mInvertSurfaceMask, mRelaxDisorientedTriangles;
255 std::unique_ptr<uint32_t[]> mQuantizedSeamPoints;
256 std::vector<uint8_t> mPointFlags;
270 const std::vector<Vec3d>&
points,
271 const std::vector<Vec3d>& normals)
277 if (points.empty())
return avgPos;
279 for (
size_t n = 0,
N = points.size();
n <
N; ++
n) {
283 avgPos /= double(points.size());
287 double m00=0,m01=0,m02=0,
294 for (
size_t n = 0, N = points.size();
n <
N; ++
n) {
296 const Vec3d& n_ref = normals[
n];
299 m00 += n_ref[0] * n_ref[0];
300 m11 += n_ref[1] * n_ref[1];
301 m22 += n_ref[2] * n_ref[2];
303 m01 += n_ref[0] * n_ref[1];
304 m02 += n_ref[0] * n_ref[2];
305 m12 += n_ref[1] * n_ref[2];
308 rhs += n_ref * n_ref.
dot(points[
n] - avgPos);
333 Mat3d D = Mat3d::identity();
341 for (
int i = 0; i < 3; ++i ) {
342 if (
std::abs(eigenValues[i]) < tolerance) {
346 D[i][i] = 1.0 / eigenValues[i];
353 return avgPos + pseudoInv * rhs;
368 namespace volume_to_mesh_internal {
370 template<
typename ValueType>
373 FillArray(ValueType*
array,
const ValueType&
v) : mArray(array), mValue(v) { }
375 void operator()(
const tbb::blocked_range<size_t>&
range)
const {
376 const ValueType
v = mValue;
377 for (
size_t n = range.begin(),
N = range.
end();
n <
N; ++
n) {
382 ValueType *
const mArray;
383 const ValueType mValue;
387 template<
typename ValueType>
389 fillArray(ValueType*
array,
const ValueType&
val,
const size_t length)
391 const auto grainSize = std::max<size_t>(
392 length / tbb::this_task_arena::max_concurrency(), 1024);
393 const tbb::blocked_range<size_t>
range(0, length, grainSize);
394 tbb::parallel_for(range, FillArray<ValueType>(array, val), tbb::simple_partitioner());
399 enum { SIGNS = 0xFF,
EDGES = 0xE00,
INSIDE = 0x100,
400 XEDGE = 0x200, YEDGE = 0x400, ZEDGE = 0x800,
SEAM = 0x1000};
404 const bool sAdaptable[256] = {
405 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,
406 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,
407 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,
408 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,
409 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,
410 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,
411 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,
412 1,0,1,0,1,0,1,0,1,0,0,0,1,0,1,1,1,1,1,1,1,0,1,1,1,1,0,1,1,1,1,1};
416 const unsigned char sAmbiguousFace[256] = {
417 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,
418 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,
419 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,
420 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,
421 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,
422 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,
423 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,
424 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
430 const unsigned char sEdgeGroupTable[256][13] = {
431 {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},
432 {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},
433 {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},
434 {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},
435 {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},
436 {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},
437 {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},
438 {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},
439 {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},
440 {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},
441 {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},
442 {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},
443 {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},
444 {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},
445 {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},
446 {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},
447 {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},
448 {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},
449 {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},
450 {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},
451 {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},
452 {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},
453 {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},
454 {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},
455 {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},
456 {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},
457 {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},
458 {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},
459 {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},
460 {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},
461 {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},
462 {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},
463 {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},
464 {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},
465 {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},
466 {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},
467 {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},
468 {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},
469 {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},
470 {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},
471 {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},
472 {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},
473 {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},
474 {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},
475 {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},
476 {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},
477 {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},
478 {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},
479 {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},
480 {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},
481 {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},
482 {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},
483 {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},
484 {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},
485 {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},
486 {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},
487 {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},
488 {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},
489 {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},
490 {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},
491 {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},
492 {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},
493 {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},
494 {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},
495 {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},
496 {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},
497 {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},
498 {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},
499 {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},
500 {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},
501 {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},
502 {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},
503 {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},
504 {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},
505 {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},
506 {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},
507 {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},
508 {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},
509 {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},
510 {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},
511 {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},
512 {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},
513 {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},
514 {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},
515 {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},
516 {0,0,0,0,0,0,0,0,0,0,0,0,0}};
525 double epsilon = 0.001)
530 const Vec3d centroid = (p0 + p1 + p2 + p3);
531 const double d = centroid.
dot(normal) * 0.25;
535 double absDist =
std::abs(p0.dot(normal) - d);
536 if (absDist > epsilon)
return false;
538 absDist =
std::abs(p1.dot(normal) - d);
539 if (absDist > epsilon)
return false;
541 absDist =
std::abs(p2.dot(normal) - d);
542 if (absDist > epsilon)
return false;
544 absDist =
std::abs(p3.dot(normal) - d);
545 if (absDist > epsilon)
return false;
558 MASK_FIRST_10_BITS = 0x000003FF,
559 MASK_DIRTY_BIT = 0x80000000,
560 MASK_INVALID_BIT = 0x40000000
564 packPoint(
const Vec3d& v)
569 assert(!(v.x() > 1.0) && !(v.y() > 1.0) && !(v.z() > 1.0));
570 assert(!(v.x() < 0.0) && !(v.y() < 0.0) && !(v.z() < 0.0));
572 data |= (uint32_t(v.x() * 1023.0) & MASK_FIRST_10_BITS) << 20;
573 data |= (uint32_t(v.y() * 1023.0) & MASK_FIRST_10_BITS) << 10;
574 data |= (uint32_t(v.z() * 1023.0) & MASK_FIRST_10_BITS);
580 unpackPoint(uint32_t data)
583 v.z() = double(data & MASK_FIRST_10_BITS) * 0.0009775171;
585 v.y() = double(data & MASK_FIRST_10_BITS) * 0.0009775171;
587 v.x() = double(data & MASK_FIRST_10_BITS) * 0.0009775171;
597 inline bool isBoolValue() {
return false; }
600 inline bool isBoolValue<bool>() {
return true; }
605 inline bool isInsideValue(T
value, T isovalue) {
return value < isovalue; }
608 inline bool isInsideValue<bool>(
bool value,
bool ) {
return value; }
611 template<
typename AccessorT>
613 getCellVertexValues(
const AccessorT& accessor, Coord ijk,
614 math::Tuple<8, typename AccessorT::ValueType>&
values)
616 values[0] = accessor.getValue(ijk);
618 values[1] = accessor.getValue(ijk);
620 values[2] = accessor.getValue(ijk);
622 values[3] = accessor.getValue(ijk);
624 values[4] = accessor.getValue(ijk);
626 values[5] = accessor.getValue(ijk);
628 values[6] = accessor.getValue(ijk);
630 values[7] = accessor.getValue(ijk);
634 template<
typename LeafT>
636 getCellVertexValues(
const LeafT& leaf,
const Index offset,
637 math::Tuple<8, typename LeafT::ValueType>& values)
639 values[0] = leaf.getValue(offset);
640 values[3] = leaf.getValue(offset + 1);
641 values[4] = leaf.getValue(offset + LeafT::DIM);
642 values[7] = leaf.getValue(offset + LeafT::DIM + 1);
643 values[1] = leaf.getValue(offset + (LeafT::DIM * LeafT::DIM));
644 values[2] = leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + 1);
645 values[5] = leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM);
646 values[6] = leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM + 1);
650 template<
typename ValueType>
652 computeSignFlags(
const math::Tuple<8, ValueType>& values,
const ValueType iso)
655 signs |= isInsideValue(values[0], iso) ? 1u : 0u;
656 signs |= isInsideValue(values[1], iso) ? 2u : 0u;
657 signs |= isInsideValue(values[2], iso) ? 4u : 0u;
658 signs |= isInsideValue(values[3], iso) ? 8u : 0u;
659 signs |= isInsideValue(values[4], iso) ? 16u : 0u;
660 signs |= isInsideValue(values[5], iso) ? 32u : 0u;
661 signs |= isInsideValue(values[6], iso) ? 64u : 0u;
662 signs |= isInsideValue(values[7], iso) ? 128u : 0u;
663 return uint8_t(signs);
669 template<
typename AccessorT>
671 evalCellSigns(
const AccessorT& accessor,
const Coord& ijk,
typename AccessorT::ValueType iso)
675 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 1u;
677 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 2u;
679 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 4u;
681 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 8u;
682 coord[1] += 1; coord[2] = ijk[2];
683 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 16u;
685 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 32u;
687 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 64u;
689 if (isInsideValue(accessor.getValue(coord), iso)) signs |= 128u;
690 return uint8_t(signs);
696 template<
typename LeafT>
698 evalCellSigns(
const LeafT& leaf,
const Index offset,
typename LeafT::ValueType iso)
703 if (isInsideValue(leaf.getValue(offset), iso)) signs |= 1u;
706 if (isInsideValue(leaf.getValue(offset + 1), iso)) signs |= 8u;
709 if (isInsideValue(leaf.getValue(offset + LeafT::DIM), iso)) signs |= 16u;
712 if (isInsideValue(leaf.getValue(offset + LeafT::DIM + 1), iso)) signs |= 128u;
715 if (isInsideValue(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) ), iso)) signs |= 2u;
718 if (isInsideValue(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + 1), iso)) signs |= 4u;
721 if (isInsideValue(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM), iso)) signs |= 32u;
724 if (isInsideValue(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM + 1), iso)) signs |= 64u;
726 return uint8_t(signs);
732 template<
class AccessorT>
734 correctCellSigns(uint8_t& signs, uint8_t
face,
735 const AccessorT& acc, Coord ijk,
typename AccessorT::ValueType iso)
740 if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 3) signs = uint8_t(~signs);
744 if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 4) signs = uint8_t(~signs);
748 if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 1) signs = uint8_t(~signs);
752 if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 2) signs = uint8_t(~signs);
756 if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 6) signs = uint8_t(~signs);
760 if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 5) signs = uint8_t(~signs);
768 template<
class AccessorT>
770 isNonManifold(
const AccessorT& accessor,
const Coord& ijk,
771 typename AccessorT::ValueType isovalue,
const int dim)
777 p[0] = isInsideValue(accessor.getValue(coord), isovalue);
779 p[1] = isInsideValue(accessor.getValue(coord), isovalue);
781 p[2] = isInsideValue(accessor.getValue(coord), isovalue);
783 p[3] = isInsideValue(accessor.getValue(coord), isovalue);
784 coord[1] += dim; coord[2] = ijk[2];
785 p[4] = isInsideValue(accessor.getValue(coord), isovalue);
787 p[5] = isInsideValue(accessor.getValue(coord), isovalue);
789 p[6] = isInsideValue(accessor.getValue(coord), isovalue);
791 p[7] = isInsideValue(accessor.getValue(coord), isovalue);
795 if (p[0]) signs |= 1u;
796 if (p[1]) signs |= 2u;
797 if (p[2]) signs |= 4u;
798 if (p[3]) signs |= 8u;
799 if (p[4]) signs |= 16u;
800 if (p[5]) signs |= 32u;
801 if (p[6]) signs |= 64u;
802 if (p[7]) signs |= 128u;
803 if (!sAdaptable[signs])
return true;
808 int i = ijk[0], ip = ijk[0] + hDim, ipp = ijk[0] + dim;
809 int j = ijk[1], jp = ijk[1] + hDim, jpp = ijk[1] + dim;
810 int k = ijk[2], kp = ijk[2] + hDim, kpp = ijk[2] + dim;
813 coord.reset(ip, j, k);
814 m = isInsideValue(accessor.getValue(coord), isovalue);
815 if (p[0] != m && p[1] != m)
return true;
818 coord.reset(ipp, j, kp);
819 m = isInsideValue(accessor.getValue(coord), isovalue);
820 if (p[1] != m && p[2] != m)
return true;
823 coord.reset(ip, j, kpp);
824 m = isInsideValue(accessor.getValue(coord), isovalue);
825 if (p[2] != m && p[3] != m)
return true;
828 coord.reset(i, j, kp);
829 m = isInsideValue(accessor.getValue(coord), isovalue);
830 if (p[0] != m && p[3] != m)
return true;
833 coord.reset(ip, jpp, k);
834 m = isInsideValue(accessor.getValue(coord), isovalue);
835 if (p[4] != m && p[5] != m)
return true;
838 coord.reset(ipp, jpp, kp);
839 m = isInsideValue(accessor.getValue(coord), isovalue);
840 if (p[5] != m && p[6] != m)
return true;
843 coord.reset(ip, jpp, kpp);
844 m = isInsideValue(accessor.getValue(coord), isovalue);
845 if (p[6] != m && p[7] != m)
return true;
848 coord.reset(i, jpp, kp);
849 m = isInsideValue(accessor.getValue(coord), isovalue);
850 if (p[7] != m && p[4] != m)
return true;
853 coord.reset(i, jp, k);
854 m = isInsideValue(accessor.getValue(coord), isovalue);
855 if (p[0] != m && p[4] != m)
return true;
858 coord.reset(ipp, jp, k);
859 m = isInsideValue(accessor.getValue(coord), isovalue);
860 if (p[1] != m && p[5] != m)
return true;
863 coord.reset(ipp, jp, kpp);
864 m = isInsideValue(accessor.getValue(coord), isovalue);
865 if (p[2] != m && p[6] != m)
return true;
869 coord.reset(i, jp, kpp);
870 m = isInsideValue(accessor.getValue(coord), isovalue);
871 if (p[3] != m && p[7] != m)
return true;
877 coord.reset(ip, jp, k);
878 m = isInsideValue(accessor.getValue(coord), isovalue);
879 if (p[0] != m && p[1] != m && p[4] != m && p[5] != m)
return true;
882 coord.reset(ipp, jp, kp);
883 m = isInsideValue(accessor.getValue(coord), isovalue);
884 if (p[1] != m && p[2] != m && p[5] != m && p[6] != m)
return true;
887 coord.reset(ip, jp, kpp);
888 m = isInsideValue(accessor.getValue(coord), isovalue);
889 if (p[2] != m && p[3] != m && p[6] != m && p[7] != m)
return true;
892 coord.reset(i, jp, kp);
893 m = isInsideValue(accessor.getValue(coord), isovalue);
894 if (p[0] != m && p[3] != m && p[4] != m && p[7] != m)
return true;
897 coord.reset(ip, j, kp);
898 m = isInsideValue(accessor.getValue(coord), isovalue);
899 if (p[0] != m && p[1] != m && p[2] != m && p[3] != m)
return true;
902 coord.reset(ip, jpp, kp);
903 m = isInsideValue(accessor.getValue(coord), isovalue);
904 if (p[4] != m && p[5] != m && p[6] != m && p[7] != m)
return true;
907 coord.reset(ip, jp, kp);
908 m = isInsideValue(accessor.getValue(coord), isovalue);
909 if (p[0] != m && p[1] != m && p[2] != m && p[3] != m &&
910 p[4] != m && p[5] != m && p[6] != m && p[7] != m)
return true;
919 template <
class LeafType>
921 mergeVoxels(LeafType& leaf,
const Coord&
start,
int dim,
int regionId)
928 for (ijk[0] = start[0]; ijk[0] < end[0]; ++ijk[0]) {
929 for (ijk[1] = start[1]; ijk[1] < end[1]; ++ijk[1]) {
930 for (ijk[2] = start[2]; ijk[2] < end[2]; ++ijk[2]) {
931 leaf.setValueOnly(ijk, regionId);
940 template <
class LeafType>
942 isMergable(LeafType& leaf,
const Coord& start,
int dim,
945 if (adaptivity < 1e-6)
return false;
947 using VecT =
typename LeafType::ValueType;
948 Coord ijk, end =
start;
953 std::vector<VecT> norms;
954 for (ijk[0] = start[0]; ijk[0] < end[0]; ++ijk[0]) {
955 for (ijk[1] = start[1]; ijk[1] < end[1]; ++ijk[1]) {
956 for (ijk[2] = start[2]; ijk[2] < end[2]; ++ijk[2]) {
958 if(!leaf.isValueOn(ijk))
continue;
959 norms.push_back(leaf.getValue(ijk));
964 size_t N = norms.size();
965 for (
size_t ni = 0; ni <
N; ++ni) {
966 VecT n_i = norms[ni];
967 for (
size_t nj = 0; nj <
N; ++nj) {
968 VecT n_j = norms[nj];
969 if ((1.0 - n_i.dot(n_j)) > adaptivity)
return false;
980 inline double evalZeroCrossing(
double v0,
double v1,
double iso) {
return (iso - v0) / (v1 -
v0); }
984 template<
typename LeafT>
986 collectCornerValues(
const LeafT& leaf,
const Index offset, std::vector<double>& values)
988 values[0] = double(leaf.getValue(offset));
989 values[3] = double(leaf.getValue(offset + 1));
990 values[4] = double(leaf.getValue(offset + LeafT::DIM));
991 values[7] = double(leaf.getValue(offset + LeafT::DIM + 1));
992 values[1] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM)));
993 values[2] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + 1));
994 values[5] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM));
995 values[6] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM + 1));
1000 template<
typename AccessorT>
1002 collectCornerValues(
const AccessorT& acc,
const Coord& ijk, std::vector<double>& values)
1005 values[0] = double(acc.getValue(coord));
1008 values[1] = double(acc.getValue(coord));
1011 values[2] = double(acc.getValue(coord));
1014 values[3] = double(acc.getValue(coord));
1016 coord[1] += 1; coord[2] = ijk[2];
1017 values[4] = double(acc.getValue(coord));
1020 values[5] = double(acc.getValue(coord));
1023 values[6] = double(acc.getValue(coord));
1026 values[7] = double(acc.getValue(coord));
1032 computePoint(
const std::vector<double>& values,
unsigned char signs,
1033 unsigned char edgeGroup,
double iso)
1035 Vec3d avg(0.0, 0.0, 0.0);
1038 if (sEdgeGroupTable[signs][1] == edgeGroup) {
1039 avg[0] += evalZeroCrossing(values[0], values[1], iso);
1043 if (sEdgeGroupTable[signs][2] == edgeGroup) {
1045 avg[2] += evalZeroCrossing(values[1], values[2], iso);
1049 if (sEdgeGroupTable[signs][3] == edgeGroup) {
1050 avg[0] += evalZeroCrossing(values[3], values[2], iso);
1055 if (sEdgeGroupTable[signs][4] == edgeGroup) {
1056 avg[2] += evalZeroCrossing(values[0], values[3], iso);
1060 if (sEdgeGroupTable[signs][5] == edgeGroup) {
1061 avg[0] += evalZeroCrossing(values[4], values[5], iso);
1066 if (sEdgeGroupTable[signs][6] == edgeGroup) {
1069 avg[2] += evalZeroCrossing(values[5], values[6], iso);
1073 if (sEdgeGroupTable[signs][7] == edgeGroup) {
1074 avg[0] += evalZeroCrossing(values[7], values[6], iso);
1080 if (sEdgeGroupTable[signs][8] == edgeGroup) {
1082 avg[2] += evalZeroCrossing(values[4], values[7], iso);
1086 if (sEdgeGroupTable[signs][9] == edgeGroup) {
1087 avg[1] += evalZeroCrossing(values[0], values[4], iso);
1091 if (sEdgeGroupTable[signs][10] == edgeGroup) {
1093 avg[1] += evalZeroCrossing(values[1], values[5], iso);
1097 if (sEdgeGroupTable[signs][11] == edgeGroup) {
1099 avg[1] += evalZeroCrossing(values[2], values[6], iso);
1104 if (sEdgeGroupTable[signs][12] == edgeGroup) {
1105 avg[1] += evalZeroCrossing(values[3], values[7], iso);
1111 double w = 1.0 / double(samples);
1124 computeMaskedPoint(
Vec3d& avg,
const std::vector<double>& values,
unsigned char signs,
1125 unsigned char signsMask,
unsigned char edgeGroup,
double iso)
1127 avg =
Vec3d(0.0, 0.0, 0.0);
1130 if (sEdgeGroupTable[signs][1] == edgeGroup
1131 && sEdgeGroupTable[signsMask][1] == 0) {
1132 avg[0] += evalZeroCrossing(values[0], values[1], iso);
1136 if (sEdgeGroupTable[signs][2] == edgeGroup
1137 && sEdgeGroupTable[signsMask][2] == 0) {
1139 avg[2] += evalZeroCrossing(values[1], values[2], iso);
1143 if (sEdgeGroupTable[signs][3] == edgeGroup
1144 && sEdgeGroupTable[signsMask][3] == 0) {
1145 avg[0] += evalZeroCrossing(values[3], values[2], iso);
1150 if (sEdgeGroupTable[signs][4] == edgeGroup
1151 && sEdgeGroupTable[signsMask][4] == 0) {
1152 avg[2] += evalZeroCrossing(values[0], values[3], iso);
1156 if (sEdgeGroupTable[signs][5] == edgeGroup
1157 && sEdgeGroupTable[signsMask][5] == 0) {
1158 avg[0] += evalZeroCrossing(values[4], values[5], iso);
1163 if (sEdgeGroupTable[signs][6] == edgeGroup
1164 && sEdgeGroupTable[signsMask][6] == 0) {
1167 avg[2] += evalZeroCrossing(values[5], values[6], iso);
1171 if (sEdgeGroupTable[signs][7] == edgeGroup
1172 && sEdgeGroupTable[signsMask][7] == 0) {
1173 avg[0] += evalZeroCrossing(values[7], values[6], iso);
1179 if (sEdgeGroupTable[signs][8] == edgeGroup
1180 && sEdgeGroupTable[signsMask][8] == 0) {
1182 avg[2] += evalZeroCrossing(values[4], values[7], iso);
1186 if (sEdgeGroupTable[signs][9] == edgeGroup
1187 && sEdgeGroupTable[signsMask][9] == 0) {
1188 avg[1] += evalZeroCrossing(values[0], values[4], iso);
1192 if (sEdgeGroupTable[signs][10] == edgeGroup
1193 && sEdgeGroupTable[signsMask][10] == 0) {
1195 avg[1] += evalZeroCrossing(values[1], values[5], iso);
1199 if (sEdgeGroupTable[signs][11] == edgeGroup
1200 && sEdgeGroupTable[signsMask][11] == 0) {
1202 avg[1] += evalZeroCrossing(values[2], values[6], iso);
1207 if (sEdgeGroupTable[signs][12] == edgeGroup
1208 && sEdgeGroupTable[signsMask][12] == 0) {
1209 avg[1] += evalZeroCrossing(values[3], values[7], iso);
1215 double w = 1.0 / double(samples);
1228 computeWeightedPoint(
const Vec3d& p,
const std::vector<double>& values,
1229 unsigned char signs,
unsigned char edgeGroup,
double iso)
1237 Vec3d avg(0.0, 0.0, 0.0);
1239 if (sEdgeGroupTable[signs][1] == edgeGroup) {
1240 avg[0] = evalZeroCrossing(values[0], values[1], iso);
1244 samples.push_back(avg);
1245 weights.push_back((avg-p).lengthSqr());
1248 if (sEdgeGroupTable[signs][2] == edgeGroup) {
1251 avg[2] = evalZeroCrossing(values[1], values[2], iso);
1253 samples.push_back(avg);
1254 weights.push_back((avg-p).lengthSqr());
1257 if (sEdgeGroupTable[signs][3] == edgeGroup) {
1258 avg[0] = evalZeroCrossing(values[3], values[2], iso);
1262 samples.push_back(avg);
1263 weights.push_back((avg-p).lengthSqr());
1266 if (sEdgeGroupTable[signs][4] == edgeGroup) {
1269 avg[2] = evalZeroCrossing(values[0], values[3], iso);
1271 samples.push_back(avg);
1272 weights.push_back((avg-p).lengthSqr());
1275 if (sEdgeGroupTable[signs][5] == edgeGroup) {
1276 avg[0] = evalZeroCrossing(values[4], values[5], iso);
1280 samples.push_back(avg);
1281 weights.push_back((avg-p).lengthSqr());
1284 if (sEdgeGroupTable[signs][6] == edgeGroup) {
1287 avg[2] = evalZeroCrossing(values[5], values[6], iso);
1289 samples.push_back(avg);
1290 weights.push_back((avg-p).lengthSqr());
1293 if (sEdgeGroupTable[signs][7] == edgeGroup) {
1294 avg[0] = evalZeroCrossing(values[7], values[6], iso);
1298 samples.push_back(avg);
1299 weights.push_back((avg-p).lengthSqr());
1302 if (sEdgeGroupTable[signs][8] == edgeGroup) {
1305 avg[2] = evalZeroCrossing(values[4], values[7], iso);
1307 samples.push_back(avg);
1308 weights.push_back((avg-p).lengthSqr());
1311 if (sEdgeGroupTable[signs][9] == edgeGroup) {
1313 avg[1] = evalZeroCrossing(values[0], values[4], iso);
1316 samples.push_back(avg);
1317 weights.push_back((avg-p).lengthSqr());
1320 if (sEdgeGroupTable[signs][10] == edgeGroup) {
1322 avg[1] = evalZeroCrossing(values[1], values[5], iso);
1325 samples.push_back(avg);
1326 weights.push_back((avg-p).lengthSqr());
1329 if (sEdgeGroupTable[signs][11] == edgeGroup) {
1331 avg[1] = evalZeroCrossing(values[2], values[6], iso);
1334 samples.push_back(avg);
1335 weights.push_back((avg-p).lengthSqr());
1338 if (sEdgeGroupTable[signs][12] == edgeGroup) {
1340 avg[1] = evalZeroCrossing(values[3], values[7], iso);
1343 samples.push_back(avg);
1344 weights.push_back((avg-p).lengthSqr());
1351 for (
size_t i = 0, I = weights.size(); i < I; ++i) {
1352 minWeight =
std::min(minWeight, weights[i]);
1353 maxWeight =
std::max(maxWeight, weights[i]);
1356 const double offset = maxWeight + minWeight * 0.1;
1357 for (
size_t i = 0, I = weights.size(); i < I; ++i) {
1358 weights[i] = offset - weights[i];
1362 double weightSum = 0.0;
1363 for (
size_t i = 0, I = weights.size(); i < I; ++i) {
1364 weightSum += weights[i];
1371 if (samples.size() > 1) {
1372 for (
size_t i = 0, I = samples.size(); i < I; ++i) {
1373 avg += samples[i] * (weights[i] / weightSum);
1376 avg = samples.front();
1386 computeCellPoints(std::vector<Vec3d>&
points,
1387 const std::vector<double>& values,
unsigned char signs,
double iso)
1389 for (
size_t n = 1, N = sEdgeGroupTable[signs][0] + 1;
n <
N; ++
n) {
1390 points.push_back(computePoint(values, signs, uint8_t(
n), iso));
1399 matchEdgeGroup(
unsigned char groupId,
unsigned char lhsSigns,
unsigned char rhsSigns)
1402 for (
size_t i = 1; i <= 12; ++i) {
1403 if (sEdgeGroupTable[lhsSigns][i] == groupId && sEdgeGroupTable[rhsSigns][i] != 0) {
1404 id = sEdgeGroupTable[rhsSigns][i];
1417 computeCellPoints(std::vector<Vec3d>& points, std::vector<bool>& weightedPointMask,
1418 const std::vector<double>& lhsValues,
const std::vector<double>& rhsValues,
1419 unsigned char lhsSigns,
unsigned char rhsSigns,
1420 double iso,
size_t pointIdx,
const uint32_t * seamPointArray)
1422 for (
size_t n = 1, N = sEdgeGroupTable[lhsSigns][0] + 1;
n <
N; ++
n) {
1424 int id = matchEdgeGroup(uint8_t(
n), lhsSigns, rhsSigns);
1428 const unsigned char e = uint8_t(
id);
1429 const uint32_t& quantizedPoint = seamPointArray[pointIdx + (
id - 1)];
1431 if ((quantizedPoint & MASK_DIRTY_BIT) && !(quantizedPoint & MASK_INVALID_BIT)) {
1432 Vec3d p = unpackPoint(quantizedPoint);
1433 points.push_back(computeWeightedPoint(p, rhsValues, rhsSigns, e, iso));
1434 weightedPointMask.push_back(
true);
1436 points.push_back(computePoint(rhsValues, rhsSigns, e, iso));
1437 weightedPointMask.push_back(
false);
1441 points.push_back(computePoint(lhsValues, lhsSigns, uint8_t(
n), iso));
1442 weightedPointMask.push_back(
false);
1448 template <
typename InputTreeType>
1449 struct ComputePoints
1451 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
1452 using InputValueType =
typename InputLeafNodeType::ValueType;
1455 using Int16LeafNodeType =
typename Int16TreeType::LeafNodeType;
1458 using Index32LeafNodeType =
typename Index32TreeType::LeafNodeType;
1460 ComputePoints(
Vec3s * pointArray,
1461 const InputTreeType& inputTree,
1462 const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
1463 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
1464 const std::unique_ptr<
Index32[]>& leafNodeOffsets,
1465 const math::Transform& xform,
1468 void setRefData(
const InputTreeType& refInputTree,
1469 const Index32TreeType& refPointIndexTree,
1470 const Int16TreeType& refSignFlagsTree,
1471 const uint32_t * quantizedSeamLinePoints,
1472 uint8_t * seamLinePointsFlags);
1474 void operator()(
const tbb::blocked_range<size_t>&)
const;
1477 Vec3s *
const mPoints;
1478 InputTreeType
const *
const mInputTree;
1479 Index32LeafNodeType *
const *
const mPointIndexNodes;
1480 Int16LeafNodeType
const *
const *
const mSignFlagsNodes;
1481 Index32 const *
const mNodeOffsets;
1482 math::Transform
const mTransform;
1483 double const mIsovalue;
1485 InputTreeType
const * mRefInputTree;
1486 Index32TreeType
const * mRefPointIndexTree;
1487 Int16TreeType
const * mRefSignFlagsTree;
1488 uint32_t
const * mQuantizedSeamLinePoints;
1489 uint8_t * mSeamLinePointsFlags;
1493 template <
typename InputTreeType>
1494 ComputePoints<InputTreeType>::ComputePoints(
1496 const InputTreeType& inputTree,
1497 const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
1498 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
1499 const std::unique_ptr<
Index32[]>& leafNodeOffsets,
1500 const math::Transform& xform,
1502 : mPoints(pointArray)
1503 , mInputTree(&inputTree)
1504 , mPointIndexNodes(pointIndexLeafNodes.data())
1505 , mSignFlagsNodes(signFlagsLeafNodes.data())
1506 , mNodeOffsets(leafNodeOffsets.get())
1509 , mRefInputTree(nullptr)
1510 , mRefPointIndexTree(nullptr)
1511 , mRefSignFlagsTree(nullptr)
1512 , mQuantizedSeamLinePoints(nullptr)
1513 , mSeamLinePointsFlags(nullptr)
1517 template <
typename InputTreeType>
1519 ComputePoints<InputTreeType>::setRefData(
1520 const InputTreeType& refInputTree,
1521 const Index32TreeType& refPointIndexTree,
1522 const Int16TreeType& refSignFlagsTree,
1523 const uint32_t * quantizedSeamLinePoints,
1524 uint8_t * seamLinePointsFlags)
1526 mRefInputTree = &refInputTree;
1527 mRefPointIndexTree = &refPointIndexTree;
1528 mRefSignFlagsTree = &refSignFlagsTree;
1529 mQuantizedSeamLinePoints = quantizedSeamLinePoints;
1530 mSeamLinePointsFlags = seamLinePointsFlags;
1533 template <
typename InputTreeType>
1535 ComputePoints<InputTreeType>::operator()(
const tbb::blocked_range<size_t>& range)
const
1537 using InputTreeAccessor = tree::ValueAccessor<const InputTreeType>;
1538 using Index32TreeAccessor = tree::ValueAccessor<const Index32TreeType>;
1539 using Int16TreeAccessor = tree::ValueAccessor<const Int16TreeType>;
1541 using IndexType =
typename Index32TreeType::ValueType;
1544 using IndexArrayMap = std::map<IndexType, IndexArray>;
1546 InputTreeAccessor inputAcc(*mInputTree);
1550 std::vector<Vec3d>
points(4);
1551 std::vector<bool> weightedPointMask(4);
1552 std::vector<double>
values(8), refValues(8);
1553 const double iso = mIsovalue;
1557 std::unique_ptr<InputTreeAccessor> refInputAcc;
1558 std::unique_ptr<Index32TreeAccessor> refPointIndexAcc;
1559 std::unique_ptr<Int16TreeAccessor> refSignFlagsAcc;
1561 const bool hasReferenceData = mRefInputTree && mRefPointIndexTree && mRefSignFlagsTree;
1563 if (hasReferenceData) {
1564 refInputAcc.reset(
new InputTreeAccessor(*mRefInputTree));
1565 refPointIndexAcc.reset(
new Index32TreeAccessor(*mRefPointIndexTree));
1566 refSignFlagsAcc.reset(
new Int16TreeAccessor(*mRefSignFlagsTree));
1569 for (
size_t n = range.begin(), N = range.end();
n !=
N; ++
n) {
1571 Index32LeafNodeType& pointIndexNode = *mPointIndexNodes[
n];
1572 const Coord& origin = pointIndexNode.origin();
1574 const Int16LeafNodeType& signFlagsNode = *mSignFlagsNodes[
n];
1575 const InputLeafNodeType * inputNode = inputAcc.probeConstLeaf(origin);
1578 const InputLeafNodeType * refInputNode =
nullptr;
1579 const Index32LeafNodeType * refPointIndexNode =
nullptr;
1580 const Int16LeafNodeType * refSignFlagsNode =
nullptr;
1582 if (hasReferenceData) {
1583 refInputNode = refInputAcc->probeConstLeaf(origin);
1584 refPointIndexNode = refPointIndexAcc->probeConstLeaf(origin);
1585 refSignFlagsNode = refSignFlagsAcc->probeConstLeaf(origin);
1588 IndexType pointOffset = IndexType(mNodeOffsets[
n]);
1589 IndexArrayMap regions;
1591 for (
auto it = pointIndexNode.beginValueOn(); it; ++it) {
1592 const Index offset = it.pos();
1594 const IndexType
id = it.getValue();
1597 regions[
id].push_back(offset);
1602 pointIndexNode.setValueOnly(offset, pointOffset);
1604 const Int16 flags = signFlagsNode.getValue(offset);
1605 uint8_t signs = uint8_t(SIGNS & flags);
1606 uint8_t refSigns = 0;
1608 if ((flags & SEAM) && refPointIndexNode && refSignFlagsNode) {
1609 if (refSignFlagsNode->isValueOn(offset)) {
1610 refSigns = uint8_t(SIGNS & refSignFlagsNode->getValue(offset));
1614 ijk = Index32LeafNodeType::offsetToLocalCoord(offset);
1616 const bool inclusiveCell = inputNode &&
1617 ijk[0] <
int(Index32LeafNodeType::DIM - 1) &&
1618 ijk[1] <
int(Index32LeafNodeType::DIM - 1) &&
1619 ijk[2] <
int(Index32LeafNodeType::DIM - 1);
1623 if (inclusiveCell) collectCornerValues(*inputNode, offset, values);
1624 else collectCornerValues(inputAcc, ijk, values);
1627 weightedPointMask.clear();
1629 if (refSigns == 0) {
1631 computeCellPoints(points, values, signs, iso);
1634 if (inclusiveCell && refInputNode) {
1635 collectCornerValues(*refInputNode, offset, refValues);
1637 collectCornerValues(*refInputAcc, ijk, refValues);
1639 computeCellPoints(points, weightedPointMask, values, refValues, signs, refSigns,
1640 iso, refPointIndexNode->getValue(offset), mQuantizedSeamLinePoints);
1643 xyz[0] = double(ijk[0]);
1644 xyz[1] = double(ijk[1]);
1645 xyz[2] = double(ijk[2]);
1647 for (
size_t i = 0, I = points.size(); i < I; ++i) {
1649 Vec3d& point = points[i];
1652 if (!std::isfinite(point[0]) ||
1653 !std::isfinite(point[1]) ||
1654 !std::isfinite(point[2]))
1657 "VolumeToMesh encountered NaNs or infs in the input VDB!"
1658 " Hint: Check the input and consider using the \"Diagnostics\" tool "
1659 "to detect and resolve the NaNs.");
1663 point = mTransform.indexToWorld(point);
1665 Vec3s& pos = mPoints[pointOffset];
1666 pos[0] = float(point[0]);
1667 pos[1] = float(point[1]);
1668 pos[2] = float(point[2]);
1670 if (mSeamLinePointsFlags && !weightedPointMask.empty() && weightedPointMask[i]) {
1671 mSeamLinePointsFlags[pointOffset] = uint8_t(1);
1679 for (
typename IndexArrayMap::iterator it = regions.begin(); it != regions.end(); ++it) {
1681 Vec3d avg(0.0), point(0.0);
1685 for (
size_t i = 0, I = voxels.size(); i < I; ++i) {
1687 const Index offset = voxels[i];
1688 ijk = Index32LeafNodeType::offsetToLocalCoord(offset);
1690 const bool inclusiveCell = inputNode &&
1691 ijk[0] <
int(Index32LeafNodeType::DIM - 1) &&
1692 ijk[1] <
int(Index32LeafNodeType::DIM - 1) &&
1693 ijk[2] <
int(Index32LeafNodeType::DIM - 1);
1697 pointIndexNode.setValueOnly(offset, pointOffset);
1699 uint8_t signs = uint8_t(SIGNS & signFlagsNode.getValue(offset));
1701 if (inclusiveCell) collectCornerValues(*inputNode, offset, values);
1702 else collectCornerValues(inputAcc, ijk, values);
1705 computeCellPoints(points, values, signs, iso);
1707 avg[0] += double(ijk[0]) + points[0][0];
1708 avg[1] += double(ijk[1]) + points[0][1];
1709 avg[2] += double(ijk[2]) + points[0][2];
1715 double w = 1.0 / double(count);
1721 avg = mTransform.indexToWorld(avg);
1723 Vec3s& pos = mPoints[pointOffset];
1724 pos[0] = float(avg[0]);
1725 pos[1] = float(avg[1]);
1726 pos[2] = float(avg[2]);
1737 template <
typename InputTreeType>
1738 struct SeamLineWeights
1740 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
1741 using InputValueType =
typename InputLeafNodeType::ValueType;
1744 using Int16LeafNodeType =
typename Int16TreeType::LeafNodeType;
1747 using Index32LeafNodeType =
typename Index32TreeType::LeafNodeType;
1749 SeamLineWeights(
const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
1750 const InputTreeType& inputTree,
1751 const Index32TreeType& refPointIndexTree,
1752 const Int16TreeType& refSignFlagsTree,
1753 uint32_t * quantizedPoints,
1755 : mSignFlagsNodes(signFlagsLeafNodes.data())
1756 , mInputTree(&inputTree)
1757 , mRefPointIndexTree(&refPointIndexTree)
1758 , mRefSignFlagsTree(&refSignFlagsTree)
1759 , mQuantizedPoints(quantizedPoints)
1764 void operator()(
const tbb::blocked_range<size_t>& range)
const
1766 tree::ValueAccessor<const InputTreeType> inputTreeAcc(*mInputTree);
1767 tree::ValueAccessor<const Index32TreeType> pointIndexTreeAcc(*mRefPointIndexTree);
1768 tree::ValueAccessor<const Int16TreeType> signFlagsTreeAcc(*mRefSignFlagsTree);
1770 std::vector<double>
values(8);
1771 const double iso = double(mIsovalue);
1775 for (
size_t n = range.begin(), N = range.end();
n !=
N; ++
n) {
1777 const Int16LeafNodeType& signFlagsNode = *mSignFlagsNodes[
n];
1778 const Coord& origin = signFlagsNode.origin();
1780 const Int16LeafNodeType * refSignNode = signFlagsTreeAcc.probeConstLeaf(origin);
1781 if (!refSignNode)
continue;
1783 const Index32LeafNodeType* refPointIndexNode =
1784 pointIndexTreeAcc.probeConstLeaf(origin);
1785 if (!refPointIndexNode)
continue;
1787 const InputLeafNodeType * inputNode = inputTreeAcc.probeConstLeaf(origin);
1789 for (
typename Int16LeafNodeType::ValueOnCIter it = signFlagsNode.cbeginValueOn();
1792 const Index offset = it.pos();
1794 ijk = Index32LeafNodeType::offsetToLocalCoord(offset);
1796 const bool inclusiveCell = inputNode &&
1797 ijk[0] <
int(Index32LeafNodeType::DIM - 1) &&
1798 ijk[1] <
int(Index32LeafNodeType::DIM - 1) &&
1799 ijk[2] <
int(Index32LeafNodeType::DIM - 1);
1803 if ((it.getValue() & SEAM) && refSignNode->isValueOn(offset)) {
1805 uint8_t lhsSigns = uint8_t(SIGNS & it.getValue());
1806 uint8_t rhsSigns = uint8_t(SIGNS & refSignNode->getValue(offset));
1809 if (inclusiveCell) {
1810 collectCornerValues(*inputNode, offset, values);
1812 collectCornerValues(inputTreeAcc, ijk, values);
1816 for (
unsigned i = 1, I = sEdgeGroupTable[lhsSigns][0] + 1; i < I; ++i) {
1818 int id = matchEdgeGroup(uint8_t(i), lhsSigns, rhsSigns);
1822 uint32_t& data = mQuantizedPoints[
1823 refPointIndexNode->getValue(offset) + (
id - 1)];
1825 if (!(data & MASK_DIRTY_BIT)) {
1827 int smaples = computeMaskedPoint(
1828 pos, values, lhsSigns, rhsSigns, uint8_t(i), iso);
1830 if (smaples > 0) data = packPoint(pos);
1831 else data = MASK_INVALID_BIT;
1833 data |= MASK_DIRTY_BIT;
1845 Int16LeafNodeType
const *
const *
const mSignFlagsNodes;
1846 InputTreeType
const *
const mInputTree;
1847 Index32TreeType
const *
const mRefPointIndexTree;
1848 Int16TreeType
const *
const mRefSignFlagsTree;
1849 uint32_t *
const mQuantizedPoints;
1850 InputValueType
const mIsovalue;
1854 template <
typename TreeType>
1855 struct SetSeamLineFlags
1857 using LeafNodeType =
typename TreeType::LeafNodeType;
1859 SetSeamLineFlags(
const std::vector<LeafNodeType*>& signFlagsLeafNodes,
1860 const TreeType& refSignFlagsTree)
1861 : mSignFlagsNodes(signFlagsLeafNodes.data())
1862 , mRefSignFlagsTree(&refSignFlagsTree)
1866 void operator()(
const tbb::blocked_range<size_t>& range)
const
1868 tree::ValueAccessor<const TreeType> refSignFlagsTreeAcc(*mRefSignFlagsTree);
1870 for (
size_t n = range.begin(), N = range.end();
n !=
N; ++
n) {
1872 LeafNodeType& signFlagsNode = *mSignFlagsNodes[
n];
1873 const Coord& origin = signFlagsNode.origin();
1875 const LeafNodeType * refSignNode = refSignFlagsTreeAcc.probeConstLeaf(origin);
1876 if (!refSignNode)
continue;
1878 for (
auto it = signFlagsNode.cbeginValueOn(); it; ++it) {
1879 const Index offset = it.pos();
1881 uint8_t rhsSigns = uint8_t(refSignNode->getValue(offset) & SIGNS);
1883 if (sEdgeGroupTable[rhsSigns][0] > 0) {
1885 const typename LeafNodeType::ValueType
value = it.getValue();
1886 uint8_t lhsSigns = uint8_t(value & SIGNS);
1888 if (rhsSigns != lhsSigns) {
1889 signFlagsNode.setValueOnly(offset, value | SEAM);
1899 LeafNodeType *
const *
const mSignFlagsNodes;
1900 TreeType
const *
const mRefSignFlagsTree;
1904 template <
typename BoolTreeType,
typename SignDataType>
1905 struct TransferSeamLineFlags
1907 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
1910 using SignDataLeafNodeType =
typename SignDataTreeType::LeafNodeType;
1912 TransferSeamLineFlags(
const std::vector<SignDataLeafNodeType*>& signFlagsLeafNodes,
1913 const BoolTreeType& maskTree)
1914 : mSignFlagsNodes(signFlagsLeafNodes.data())
1915 , mMaskTree(&maskTree)
1919 void operator()(
const tbb::blocked_range<size_t>& range)
const
1921 tree::ValueAccessor<const BoolTreeType> maskAcc(*mMaskTree);
1923 for (
size_t n = range.begin(), N = range.end();
n !=
N; ++
n) {
1925 SignDataLeafNodeType& signFlagsNode = *mSignFlagsNodes[
n];
1926 const Coord& origin = signFlagsNode.origin();
1928 const BoolLeafNodeType * maskNode = maskAcc.probeConstLeaf(origin);
1929 if (!maskNode)
continue;
1931 using ValueOnCIter =
typename SignDataLeafNodeType::ValueOnCIter;
1933 for (ValueOnCIter it = signFlagsNode.cbeginValueOn(); it; ++it) {
1934 const Index offset = it.pos();
1936 if (maskNode->isValueOn(offset)) {
1937 signFlagsNode.setValueOnly(offset, it.getValue() | SEAM);
1946 SignDataLeafNodeType *
const *
const mSignFlagsNodes;
1947 BoolTreeType
const *
const mMaskTree;
1951 template <
typename TreeType>
1952 struct MaskSeamLineVoxels
1954 using LeafNodeType =
typename TreeType::LeafNodeType;
1957 MaskSeamLineVoxels(
const std::vector<LeafNodeType*>& signFlagsLeafNodes,
1958 const TreeType& signFlagsTree,
1960 : mSignFlagsNodes(signFlagsLeafNodes.data())
1961 , mSignFlagsTree(&signFlagsTree)
1967 MaskSeamLineVoxels(MaskSeamLineVoxels& rhs,
tbb::split)
1968 : mSignFlagsNodes(rhs.mSignFlagsNodes)
1969 , mSignFlagsTree(rhs.mSignFlagsTree)
1975 void join(MaskSeamLineVoxels& rhs) { mMask->merge(*rhs.mMask); }
1977 void operator()(
const tbb::blocked_range<size_t>& range)
1979 using ValueOnCIter =
typename LeafNodeType::ValueOnCIter;
1980 using ValueType =
typename LeafNodeType::ValueType;
1982 tree::ValueAccessor<const TreeType> signFlagsAcc(*mSignFlagsTree);
1983 tree::ValueAccessor<BoolTreeType> maskAcc(*mMask);
1986 for (
size_t n = range.begin(), N = range.end();
n !=
N; ++
n) {
1988 LeafNodeType& signFlagsNode = *mSignFlagsNodes[
n];
1991 for (ValueOnCIter it = signFlagsNode.cbeginValueOn(); it; ++it) {
1993 const ValueType flags = it.getValue();
1995 if (!(flags & SEAM) && (flags & EDGES)) {
1997 ijk = it.getCoord();
1999 bool isSeamLineVoxel =
false;
2001 if (flags & XEDGE) {
2003 isSeamLineVoxel = (signFlagsAcc.getValue(ijk) & SEAM);
2005 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
2007 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
2011 if (!isSeamLineVoxel && flags & YEDGE) {
2013 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
2015 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
2017 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
2021 if (!isSeamLineVoxel && flags & ZEDGE) {
2023 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
2025 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
2027 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
2031 if (isSeamLineVoxel) {
2032 maskAcc.setValue(it.getCoord(),
true);
2041 LeafNodeType *
const *
const mSignFlagsNodes;
2042 TreeType
const *
const mSignFlagsTree;
2043 BoolTreeType mTempMask;
2044 BoolTreeType *
const mMask;
2048 template<
typename SignDataTreeType>
2050 markSeamLineData(SignDataTreeType& signFlagsTree,
const SignDataTreeType& refSignFlagsTree)
2052 using SignDataType =
typename SignDataTreeType::ValueType;
2053 using SignDataLeafNodeType =
typename SignDataTreeType::LeafNodeType;
2056 std::vector<SignDataLeafNodeType*> signFlagsLeafNodes;
2057 signFlagsTree.getNodes(signFlagsLeafNodes);
2059 const tbb::blocked_range<size_t> nodeRange(0, signFlagsLeafNodes.size());
2062 SetSeamLineFlags<SignDataTreeType>(signFlagsLeafNodes, refSignFlagsTree));
2064 BoolTreeType seamLineMaskTree(
false);
2066 MaskSeamLineVoxels<SignDataTreeType>
2067 maskSeamLine(signFlagsLeafNodes, signFlagsTree, seamLineMaskTree);
2069 tbb::parallel_reduce(nodeRange, maskSeamLine);
2072 TransferSeamLineFlags<BoolTreeType, SignDataType>(signFlagsLeafNodes, seamLineMaskTree));
2079 template <
typename InputGr
idType>
2080 struct MergeVoxelRegions
2082 using InputTreeType =
typename InputGridType::TreeType;
2083 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
2084 using InputValueType =
typename InputLeafNodeType::ValueType;
2087 using FloatLeafNodeType =
typename FloatTreeType::LeafNodeType;
2088 using FloatGridType = Grid<FloatTreeType>;
2091 using Int16LeafNodeType =
typename Int16TreeType::LeafNodeType;
2094 using Index32LeafNodeType =
typename Index32TreeType::LeafNodeType;
2096 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
2097 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
2099 MergeVoxelRegions(
const InputGridType& inputGrid,
2100 const Index32TreeType& pointIndexTree,
2101 const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
2102 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
2105 bool invertSurfaceOrientation);
2107 void setSpatialAdaptivity(
const FloatGridType& grid)
2109 mSpatialAdaptivityTree = &grid.tree();
2110 mSpatialAdaptivityTransform = &grid.transform();
2113 void setAdaptivityMask(
const BoolTreeType&
mask)
2118 void setRefSignFlagsData(
const Int16TreeType& signFlagsData,
float internalAdaptivity)
2120 mRefSignFlagsTree = &signFlagsData;
2121 mInternalAdaptivity = internalAdaptivity;
2124 void operator()(
const tbb::blocked_range<size_t>&)
const;
2127 InputTreeType
const *
const mInputTree;
2128 math::Transform
const *
const mInputTransform;
2130 Index32TreeType
const *
const mPointIndexTree;
2131 Index32LeafNodeType *
const *
const mPointIndexNodes;
2132 Int16LeafNodeType
const *
const *
const mSignFlagsNodes;
2134 InputValueType mIsovalue;
2135 float mSurfaceAdaptivity, mInternalAdaptivity;
2136 bool mInvertSurfaceOrientation;
2138 FloatTreeType
const * mSpatialAdaptivityTree;
2139 BoolTreeType
const * mMaskTree;
2140 Int16TreeType
const * mRefSignFlagsTree;
2141 math::Transform
const * mSpatialAdaptivityTransform;
2145 template <
typename InputGr
idType>
2146 MergeVoxelRegions<InputGridType>::MergeVoxelRegions(
2147 const InputGridType& inputGrid,
2148 const Index32TreeType& pointIndexTree,
2149 const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
2150 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
2153 bool invertSurfaceOrientation)
2154 : mInputTree(&inputGrid.tree())
2155 , mInputTransform(&inputGrid.
transform())
2156 , mPointIndexTree(&pointIndexTree)
2157 , mPointIndexNodes(pointIndexLeafNodes.data())
2158 , mSignFlagsNodes(signFlagsLeafNodes.data())
2160 , mSurfaceAdaptivity(adaptivity)
2161 , mInternalAdaptivity(adaptivity)
2162 , mInvertSurfaceOrientation(invertSurfaceOrientation)
2163 , mSpatialAdaptivityTree(nullptr)
2164 , mMaskTree(nullptr)
2165 , mRefSignFlagsTree(nullptr)
2166 , mSpatialAdaptivityTransform(nullptr)
2171 template <
typename InputGr
idType>
2173 MergeVoxelRegions<InputGridType>::operator()(
const tbb::blocked_range<size_t>& range)
const
2175 using Vec3sType = math::Vec3<float>;
2178 using InputTreeAccessor = tree::ValueAccessor<const InputTreeType>;
2179 using FloatTreeAccessor = tree::ValueAccessor<const FloatTreeType>;
2180 using Index32TreeAccessor = tree::ValueAccessor<const Index32TreeType>;
2181 using Int16TreeAccessor = tree::ValueAccessor<const Int16TreeType>;
2182 using BoolTreeAccessor = tree::ValueAccessor<const BoolTreeType>;
2184 std::unique_ptr<FloatTreeAccessor> spatialAdaptivityAcc;
2185 if (mSpatialAdaptivityTree && mSpatialAdaptivityTransform) {
2186 spatialAdaptivityAcc.reset(
new FloatTreeAccessor(*mSpatialAdaptivityTree));
2189 std::unique_ptr<BoolTreeAccessor> maskAcc;
2191 maskAcc.reset(
new BoolTreeAccessor(*mMaskTree));
2194 std::unique_ptr<Int16TreeAccessor> refSignFlagsAcc;
2195 if (mRefSignFlagsTree) {
2196 refSignFlagsAcc.reset(
new Int16TreeAccessor(*mRefSignFlagsTree));
2199 InputTreeAccessor inputAcc(*mInputTree);
2200 Index32TreeAccessor pointIndexAcc(*mPointIndexTree);
2202 BoolLeafNodeType
mask;
2204 const bool invertGradientDir = mInvertSurfaceOrientation || isBoolValue<InputValueType>();
2205 std::unique_ptr<Vec3sLeafNodeType> gradientNode;
2208 const int LeafDim = InputLeafNodeType::DIM;
2210 for (
size_t n = range.begin(), N = range.end();
n !=
N; ++
n) {
2212 mask.setValuesOff();
2214 const Int16LeafNodeType& signFlagsNode = *mSignFlagsNodes[
n];
2215 Index32LeafNodeType& pointIndexNode = *mPointIndexNodes[
n];
2217 const Coord& origin = pointIndexNode.origin();
2219 end[0] = origin[0] + LeafDim;
2220 end[1] = origin[1] + LeafDim;
2221 end[2] = origin[2] + LeafDim;
2225 const BoolLeafNodeType* maskLeaf = maskAcc->probeConstLeaf(origin);
2226 if (maskLeaf !=
nullptr) {
2227 for (
typename BoolLeafNodeType::ValueOnCIter it = maskLeaf->cbeginValueOn();
2230 mask.setActiveState(it.getCoord() & ~1u,
true);
2235 float adaptivity = (refSignFlagsAcc && !refSignFlagsAcc->probeConstLeaf(origin)) ?
2236 mInternalAdaptivity : mSurfaceAdaptivity;
2238 bool useGradients = adaptivity < 1.0f;
2241 FloatLeafNodeType adaptivityLeaf(origin, adaptivity);
2243 if (spatialAdaptivityAcc) {
2244 useGradients =
false;
2245 for (
Index offset = 0; offset < FloatLeafNodeType::NUM_VALUES; ++
offset) {
2246 ijk = adaptivityLeaf.offsetToGlobalCoord(offset);
2247 ijk = mSpatialAdaptivityTransform->worldToIndexCellCentered(
2248 mInputTransform->indexToWorld(ijk));
2249 float weight = spatialAdaptivityAcc->getValue(ijk);
2250 float adaptivityValue = weight * adaptivity;
2251 if (adaptivityValue < 1.0
f) useGradients =
true;
2252 adaptivityLeaf.setValueOnly(offset, adaptivityValue);
2257 for (
auto it = signFlagsNode.cbeginValueOn(); it; ++it) {
2258 const Int16 flags = it.getValue();
2259 const unsigned char signs =
static_cast<unsigned char>(SIGNS &
int(flags));
2261 if ((flags & SEAM) || !sAdaptable[signs] || sEdgeGroupTable[signs][0] > 1) {
2263 mask.setActiveState(it.getCoord() & ~1u,
true);
2265 }
else if (flags & EDGES) {
2267 bool maskRegion =
false;
2269 ijk = it.getCoord();
2270 if (!pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2272 if (!maskRegion && flags & XEDGE) {
2274 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2276 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2278 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2282 if (!maskRegion && flags & YEDGE) {
2284 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2286 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2288 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2292 if (!maskRegion && flags & ZEDGE) {
2294 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2296 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2298 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2303 mask.setActiveState(it.getCoord() & ~1u,
true);
2310 for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
2311 for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
2312 for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
2313 if (!mask.isValueOn(ijk) && isNonManifold(inputAcc, ijk, mIsovalue, dim)) {
2314 mask.setActiveState(ijk,
true);
2325 gradientNode->setValuesOff();
2327 gradientNode.reset(
new Vec3sLeafNodeType());
2330 for (
auto it = signFlagsNode.cbeginValueOn(); it; ++it) {
2331 ijk = it.getCoord();
2332 if (!mask.isValueOn(ijk & ~1u)) {
2336 if (invertGradientDir) {
2340 gradientNode->setValueOn(it.pos(), dir);
2347 for ( ; dim <= LeafDim; dim = dim << 1) {
2348 const unsigned coordMask = ~((dim << 1) - 1);
2349 for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
2350 for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
2351 for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
2353 adaptivity = adaptivityLeaf.getValue(ijk);
2355 if (mask.isValueOn(ijk)
2356 || isNonManifold(inputAcc, ijk, mIsovalue, dim)
2357 || (useGradients && !isMergable(*gradientNode, ijk, dim, adaptivity)))
2359 mask.setActiveState(ijk & coordMask,
true);
2361 mergeVoxels(pointIndexNode, ijk, dim, regionId++);
2375 struct UniformPrimBuilder
2377 UniformPrimBuilder(): mIdx(0), mPolygonPool(nullptr) {}
2379 void init(
const size_t upperBound, PolygonPool& quadPool)
2381 mPolygonPool = &quadPool;
2382 mPolygonPool->resetQuads(upperBound);
2386 template<
typename IndexType>
2387 void addPrim(
const math::Vec4<IndexType>& verts,
bool reverse,
char flags = 0)
2390 mPolygonPool->quad(mIdx) = verts;
2392 Vec4I& quad = mPolygonPool->quad(mIdx);
2398 mPolygonPool->quadFlags(mIdx) =
flags;
2404 mPolygonPool->trimQuads(mIdx);
2409 PolygonPool* mPolygonPool;
2414 struct AdaptivePrimBuilder
2416 AdaptivePrimBuilder() : mQuadIdx(0), mTriangleIdx(0), mPolygonPool(nullptr) {}
2418 void init(
const size_t upperBound, PolygonPool& polygonPool)
2420 mPolygonPool = &polygonPool;
2421 mPolygonPool->resetQuads(upperBound);
2422 mPolygonPool->resetTriangles(upperBound);
2428 template<
typename IndexType>
2429 void addPrim(
const math::Vec4<IndexType>& verts,
bool reverse,
char flags = 0)
2431 if (verts[0] != verts[1] && verts[0] != verts[2] && verts[0] != verts[3]
2432 && verts[1] != verts[2] && verts[1] != verts[3] && verts[2] != verts[3]) {
2433 mPolygonPool->quadFlags(mQuadIdx) =
flags;
2434 addQuad(verts, reverse);
2436 verts[0] == verts[3] &&
2437 verts[1] != verts[2] &&
2438 verts[1] != verts[0] &&
2439 verts[2] != verts[0]) {
2440 mPolygonPool->triangleFlags(mTriangleIdx) =
flags;
2441 addTriangle(verts[0], verts[1], verts[2], reverse);
2443 verts[1] == verts[2] &&
2444 verts[0] != verts[3] &&
2445 verts[0] != verts[1] &&
2446 verts[3] != verts[1]) {
2447 mPolygonPool->triangleFlags(mTriangleIdx) =
flags;
2448 addTriangle(verts[0], verts[1], verts[3], reverse);
2450 verts[0] == verts[1] &&
2451 verts[2] != verts[3] &&
2452 verts[2] != verts[0] &&
2453 verts[3] != verts[0]) {
2454 mPolygonPool->triangleFlags(mTriangleIdx) =
flags;
2455 addTriangle(verts[0], verts[2], verts[3], reverse);
2457 verts[2] == verts[3] &&
2458 verts[0] != verts[1] &&
2459 verts[0] != verts[2] &&
2460 verts[1] != verts[2]) {
2461 mPolygonPool->triangleFlags(mTriangleIdx) =
flags;
2462 addTriangle(verts[0], verts[1], verts[2], reverse);
2469 mPolygonPool->trimQuads(mQuadIdx,
true);
2470 mPolygonPool->trimTrinagles(mTriangleIdx,
true);
2475 template<
typename IndexType>
2476 void addQuad(
const math::Vec4<IndexType>& verts,
bool reverse)
2479 mPolygonPool->quad(mQuadIdx) = verts;
2481 Vec4I& quad = mPolygonPool->quad(mQuadIdx);
2490 void addTriangle(
unsigned v0,
unsigned v1,
unsigned v2,
bool reverse)
2492 Vec3I& prim = mPolygonPool->triangle(mTriangleIdx);
2506 size_t mQuadIdx, mTriangleIdx;
2507 PolygonPool *mPolygonPool;
2511 template<
typename SignAccT,
typename IdxAccT,
typename PrimBuilder>
2514 bool invertSurfaceOrientation,
2519 const SignAccT& signAcc,
2520 const IdxAccT& idxAcc,
2521 PrimBuilder& mesher)
2523 using IndexType =
typename IdxAccT::ValueType;
2526 const bool isActive = idxAcc.probeValue(ijk, v0);
2533 bool isInside = flags &
INSIDE;
2535 isInside = invertSurfaceOrientation ? !isInside : isInside;
2538 math::Vec4<IndexType> quad(0,0,0,0);
2540 if (flags & XEDGE) {
2542 quad[0] = v0 + offsets[0];
2546 bool activeValues = idxAcc.probeValue(coord, quad[1]);
2547 uint8_t cell = uint8_t(SIGNS & signAcc.getValue(coord));
2548 quad[1] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][5] - 1 : 0;
2552 activeValues = activeValues && idxAcc.probeValue(coord, quad[2]);
2553 cell = uint8_t(SIGNS & signAcc.getValue(coord));
2554 quad[2] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][7] - 1 : 0;
2558 activeValues = activeValues && idxAcc.probeValue(coord, quad[3]);
2559 cell = uint8_t(SIGNS & signAcc.getValue(coord));
2560 quad[3] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][3] - 1 : 0;
2563 mesher.addPrim(quad, isInside, tag[
bool(refFlags & XEDGE)]);
2570 if (flags & YEDGE) {
2572 quad[0] = v0 + offsets[1];
2576 bool activeValues = idxAcc.probeValue(coord, quad[1]);
2577 uint8_t cell = uint8_t(SIGNS & signAcc.getValue(coord));
2578 quad[1] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][12] - 1 : 0;
2582 activeValues = activeValues && idxAcc.probeValue(coord, quad[2]);
2583 cell = uint8_t(SIGNS & signAcc.getValue(coord));
2584 quad[2] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][11] - 1 : 0;
2588 activeValues = activeValues && idxAcc.probeValue(coord, quad[3]);
2589 cell = uint8_t(SIGNS & signAcc.getValue(coord));
2590 quad[3] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][10] - 1 : 0;
2593 mesher.addPrim(quad, isInside, tag[
bool(refFlags & YEDGE)]);
2600 if (flags & ZEDGE) {
2602 quad[0] = v0 + offsets[2];
2606 bool activeValues = idxAcc.probeValue(coord, quad[1]);
2607 uint8_t cell = uint8_t(SIGNS & signAcc.getValue(coord));
2608 quad[1] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][8] - 1 : 0;
2612 activeValues = activeValues && idxAcc.probeValue(coord, quad[2]);
2613 cell = uint8_t(SIGNS & signAcc.getValue(coord));
2614 quad[2] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][6] - 1 : 0;
2618 activeValues = activeValues && idxAcc.probeValue(coord, quad[3]);
2619 cell = uint8_t(SIGNS & signAcc.getValue(coord));
2620 quad[3] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][2] - 1 : 0;
2623 mesher.addPrim(quad, !isInside, tag[
bool(refFlags & ZEDGE)]);
2632 template<
typename InputTreeType>
2633 struct MaskTileBorders
2635 using InputValueType =
typename InputTreeType::ValueType;
2636 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
2639 MaskTileBorders(
const InputTreeType& inputTree, InputValueType iso,
2640 BoolTreeType& mask,
const Vec4i* tileArray)
2641 : mInputTree(&inputTree)
2645 , mTileArray(tileArray)
2649 MaskTileBorders(MaskTileBorders& rhs,
tbb::split)
2650 : mInputTree(rhs.mInputTree)
2651 , mIsovalue(rhs.mIsovalue)
2654 , mTileArray(rhs.mTileArray)
2658 void join(MaskTileBorders& rhs) { mMask->merge(*rhs.mMask); }
2660 void operator()(
const tbb::blocked_range<size_t>&);
2663 InputTreeType
const *
const mInputTree;
2664 InputValueType
const mIsovalue;
2665 BoolTreeType mTempMask;
2666 BoolTreeType *
const mMask;
2667 Vec4i const *
const mTileArray;
2671 template<
typename InputTreeType>
2673 MaskTileBorders<InputTreeType>::operator()(
const tbb::blocked_range<size_t>& range)
2675 tree::ValueAccessor<const InputTreeType> inputTreeAcc(*mInputTree);
2677 CoordBBox region, bbox;
2680 for (
size_t n = range.begin(), N = range.end();
n !=
N; ++
n) {
2682 const Vec4i& tile = mTileArray[
n];
2684 bbox.min()[0] = tile[0];
2685 bbox.min()[1] = tile[1];
2686 bbox.min()[2] = tile[2];
2687 bbox.max() = bbox.min();
2688 bbox.max().offset(tile[3]);
2690 InputValueType
value = mInputTree->background();
2692 const bool isInside = isInsideValue(inputTreeAcc.getValue(bbox.min()), mIsovalue);
2693 const int valueDepth = inputTreeAcc.getValueDepth(bbox.min());
2701 bool processRegion =
true;
2702 if (valueDepth >= inputTreeAcc.getValueDepth(nijk)) {
2703 processRegion = isInside != isInsideValue(inputTreeAcc.getValue(nijk), mIsovalue);
2706 if (processRegion) {
2709 region.min()[0] = region.max()[0] = ijk[0];
2710 mMask->fill(region,
false);
2717 processRegion =
true;
2718 if (valueDepth >= inputTreeAcc.getValueDepth(ijk)) {
2719 processRegion = (!inputTreeAcc.probeValue(ijk, value)
2720 && isInside != isInsideValue(value, mIsovalue));
2723 if (processRegion) {
2726 region.min()[0] = region.max()[0] = ijk[0];
2727 mMask->fill(region,
false);
2737 processRegion =
true;
2738 if (valueDepth >= inputTreeAcc.getValueDepth(nijk)) {
2739 processRegion = isInside != isInsideValue(inputTreeAcc.getValue(nijk), mIsovalue);
2742 if (processRegion) {
2745 region.min()[1] = region.max()[1] = ijk[1];
2746 mMask->fill(region,
false);
2753 processRegion =
true;
2754 if (valueDepth >= inputTreeAcc.getValueDepth(ijk)) {
2755 processRegion = (!inputTreeAcc.probeValue(ijk, value)
2756 && isInside != isInsideValue(value, mIsovalue));
2759 if (processRegion) {
2762 region.min()[1] = region.max()[1] = ijk[1];
2763 mMask->fill(region,
false);
2773 processRegion =
true;
2774 if (valueDepth >= inputTreeAcc.getValueDepth(nijk)) {
2775 processRegion = isInside != isInsideValue(inputTreeAcc.getValue(nijk), mIsovalue);
2778 if (processRegion) {
2781 region.min()[2] = region.max()[2] = ijk[2];
2782 mMask->fill(region,
false);
2788 processRegion =
true;
2789 if (valueDepth >= inputTreeAcc.getValueDepth(ijk)) {
2790 processRegion = (!inputTreeAcc.probeValue(ijk, value)
2791 && isInside != isInsideValue(value, mIsovalue));
2794 if (processRegion) {
2797 region.min()[2] = region.max()[2] = ijk[2];
2798 mMask->fill(region,
false);
2804 template<
typename InputTreeType>
2806 maskActiveTileBorders(
const InputTreeType& inputTree,
typename InputTreeType::ValueType iso,
2807 typename InputTreeType::template ValueConverter<bool>::Type& mask)
2809 typename InputTreeType::ValueOnCIter tileIter(inputTree);
2810 tileIter.setMaxDepth(InputTreeType::ValueOnCIter::LEAF_DEPTH - 1);
2812 size_t tileCount = 0;
2813 for ( ; tileIter; ++tileIter) {
2817 if (tileCount > 0) {
2818 std::unique_ptr<Vec4i[]> tiles(
new Vec4i[tileCount]);
2823 tileIter = inputTree.cbeginValueOn();
2824 tileIter.setMaxDepth(InputTreeType::ValueOnCIter::LEAF_DEPTH - 1);
2826 for (; tileIter; ++tileIter) {
2827 Vec4i& tile = tiles[index++];
2828 tileIter.getBoundingBox(bbox);
2829 tile[0] = bbox.min()[0];
2830 tile[1] = bbox.min()[1];
2831 tile[2] = bbox.min()[2];
2832 tile[3] = bbox.max()[0] - bbox.min()[0];
2835 MaskTileBorders<InputTreeType> op(inputTree, iso, mask, tiles.get());
2836 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, tileCount), op);
2848 PointListCopy(
const PointList& pointsIn, std::vector<Vec3s>& pointsOut)
2849 : mPointsIn(pointsIn) , mPointsOut(pointsOut)
2853 void operator()(
const tbb::blocked_range<size_t>& range)
const
2855 for (
size_t n = range.begin();
n < range.end(); ++
n) {
2856 mPointsOut[
n] = mPointsIn[
n];
2862 std::vector<Vec3s>& mPointsOut;
2870 struct LeafNodeVoxelOffsets
2872 using IndexVector = std::vector<Index>;
2874 template<
typename LeafNodeType>
2875 void constructOffsetList();
2878 const IndexVector& core()
const {
return mCore; }
2882 const IndexVector& minX()
const {
return mMinX; }
2885 const IndexVector&
maxX()
const {
return mMaxX; }
2889 const IndexVector&
minY()
const {
return mMinY; }
2892 const IndexVector&
maxY()
const {
return mMaxY; }
2896 const IndexVector&
minZ()
const {
return mMinZ; }
2899 const IndexVector&
maxZ()
const {
return mMaxZ; }
2903 const IndexVector& internalNeighborsX()
const {
return mInternalNeighborsX; }
2906 const IndexVector& internalNeighborsY()
const {
return mInternalNeighborsY; }
2909 const IndexVector& internalNeighborsZ()
const {
return mInternalNeighborsZ; }
2913 IndexVector mCore, mMinX, mMaxX, mMinY, mMaxY, mMinZ, mMaxZ,
2914 mInternalNeighborsX, mInternalNeighborsY, mInternalNeighborsZ;
2918 template<
typename LeafNodeType>
2920 LeafNodeVoxelOffsets::constructOffsetList()
2924 mCore.reserve((LeafNodeType::DIM - 2) * (LeafNodeType::DIM - 2));
2926 for (
Index x = 1;
x < (LeafNodeType::DIM - 1); ++
x) {
2927 const Index offsetX =
x << (2 * LeafNodeType::LOG2DIM);
2928 for (
Index y = 1;
y < (LeafNodeType::DIM - 1); ++
y) {
2929 const Index offsetXY = offsetX + (
y << LeafNodeType::LOG2DIM);
2930 for (
Index z = 1;
z < (LeafNodeType::DIM - 1); ++
z) {
2931 mCore.push_back(offsetXY +
z);
2937 mInternalNeighborsX.clear();
2938 mInternalNeighborsX.reserve(
LeafNodeType::SIZE - (LeafNodeType::DIM * LeafNodeType::DIM));
2940 for (
Index x = 0;
x < (LeafNodeType::DIM - 1); ++
x) {
2941 const Index offsetX =
x << (2 * LeafNodeType::LOG2DIM);
2942 for (
Index y = 0;
y < LeafNodeType::DIM; ++
y) {
2943 const Index offsetXY = offsetX + (
y << LeafNodeType::LOG2DIM);
2944 for (
Index z = 0;
z < LeafNodeType::DIM; ++
z) {
2945 mInternalNeighborsX.push_back(offsetXY +
z);
2951 mInternalNeighborsY.clear();
2952 mInternalNeighborsY.reserve(
LeafNodeType::SIZE - (LeafNodeType::DIM * LeafNodeType::DIM));
2954 for (
Index x = 0;
x < LeafNodeType::DIM; ++
x) {
2955 const Index offsetX =
x << (2 * LeafNodeType::LOG2DIM);
2956 for (
Index y = 0;
y < (LeafNodeType::DIM - 1); ++
y) {
2957 const Index offsetXY = offsetX + (
y << LeafNodeType::LOG2DIM);
2958 for (
Index z = 0;
z < LeafNodeType::DIM; ++
z) {
2959 mInternalNeighborsY.push_back(offsetXY +
z);
2965 mInternalNeighborsZ.clear();
2966 mInternalNeighborsZ.reserve(
LeafNodeType::SIZE - (LeafNodeType::DIM * LeafNodeType::DIM));
2968 for (
Index x = 0;
x < LeafNodeType::DIM; ++
x) {
2969 const Index offsetX =
x << (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 - 1); ++
z) {
2973 mInternalNeighborsZ.push_back(offsetXY +
z);
2980 mMinX.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
2982 for (
Index y = 0;
y < LeafNodeType::DIM; ++
y) {
2983 const Index offsetXY = (
y << LeafNodeType::LOG2DIM);
2984 for (
Index z = 0;
z < LeafNodeType::DIM; ++
z) {
2985 mMinX.push_back(offsetXY +
z);
2992 mMaxX.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
2994 const Index offsetX = (LeafNodeType::DIM - 1) << (2 * LeafNodeType::LOG2DIM);
2995 for (
Index y = 0;
y < LeafNodeType::DIM; ++
y) {
2996 const Index offsetXY = offsetX + (
y << LeafNodeType::LOG2DIM);
2997 for (
Index z = 0;
z < LeafNodeType::DIM; ++
z) {
2998 mMaxX.push_back(offsetXY +
z);
3005 mMinY.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 z = 0;
z < (LeafNodeType::DIM - 1); ++
z) {
3010 mMinY.push_back(offsetX +
z);
3017 mMaxY.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
3019 const Index offsetY = (LeafNodeType::DIM - 1) << LeafNodeType::LOG2DIM;
3020 for (
Index x = 0;
x < LeafNodeType::DIM; ++
x) {
3021 const Index offsetX =
x << (2 * LeafNodeType::LOG2DIM);
3022 for (
Index z = 0;
z < (LeafNodeType::DIM - 1); ++
z) {
3023 mMaxY.push_back(offsetX + offsetY +
z);
3030 mMinZ.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
3032 for (
Index x = 0;
x < LeafNodeType::DIM; ++
x) {
3033 const Index offsetX =
x << (2 * LeafNodeType::LOG2DIM);
3034 for (
Index y = 0;
y < LeafNodeType::DIM; ++
y) {
3035 const Index offsetXY = offsetX + (
y << LeafNodeType::LOG2DIM);
3036 mMinZ.push_back(offsetXY);
3043 mMaxZ.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
3045 for (
Index x = 0;
x < LeafNodeType::DIM; ++
x) {
3046 const Index offsetX =
x << (2 * LeafNodeType::LOG2DIM);
3047 for (
Index y = 0;
y < LeafNodeType::DIM; ++
y) {
3048 const Index offsetXY = offsetX + (
y << LeafNodeType::LOG2DIM);
3049 mMaxZ.push_back(offsetXY + (LeafNodeType::DIM - 1));
3060 template<
typename AccessorT,
int _AXIS>
3061 struct VoxelEdgeAccessor {
3063 enum {
AXIS = _AXIS };
3066 VoxelEdgeAccessor(AccessorT& _acc) : acc(_acc) {}
3068 void set(Coord ijk) {
3070 acc.setActiveState(ijk);
3072 acc.setActiveState(ijk);
3074 acc.setActiveState(ijk);
3076 acc.setActiveState(ijk);
3077 }
else if (_AXIS == 1) {
3078 acc.setActiveState(ijk);
3080 acc.setActiveState(ijk);
3082 acc.setActiveState(ijk);
3084 acc.setActiveState(ijk);
3086 acc.setActiveState(ijk);
3088 acc.setActiveState(ijk);
3090 acc.setActiveState(ijk);
3092 acc.setActiveState(ijk);
3101 template<
typename VoxelEdgeAcc,
typename LeafNode>
3103 evalInternalVoxelEdges(VoxelEdgeAcc& edgeAcc,
const LeafNode& leafnode,
3104 const LeafNodeVoxelOffsets& voxels,
const typename LeafNode::ValueType iso)
3107 const std::vector<Index>* offsets = &voxels.internalNeighborsZ();
3109 if (VoxelEdgeAcc::AXIS == 0) {
3110 nvo = LeafNode::DIM * LeafNode::DIM;
3111 offsets = &voxels.internalNeighborsX();
3112 }
else if (VoxelEdgeAcc::AXIS == 1) {
3113 nvo = LeafNode::DIM;
3114 offsets = &voxels.internalNeighborsY();
3117 for (
size_t n = 0, N = offsets->size();
n <
N; ++
n) {
3118 const Index& pos = (*offsets)[
n];
3119 bool isActive = leafnode.isValueOn(pos) || leafnode.isValueOn(pos + nvo);
3120 if (isActive && (isInsideValue(leafnode.getValue(pos), iso) !=
3121 isInsideValue(leafnode.getValue(pos + nvo), iso))) {
3122 edgeAcc.set(leafnode.offsetToGlobalCoord(pos));
3131 template<
typename LeafNode,
typename TreeAcc,
typename VoxelEdgeAcc>
3133 evalExtrenalVoxelEdges(VoxelEdgeAcc& edgeAcc, TreeAcc& acc,
const LeafNode& lhsNode,
3134 const LeafNodeVoxelOffsets& voxels,
const typename LeafNode::ValueType iso)
3136 const std::vector<Index>* lhsOffsets = &voxels.maxX();
3137 const std::vector<Index>* rhsOffsets = &voxels.minX();
3138 Coord ijk = lhsNode.origin();
3140 if (VoxelEdgeAcc::AXIS == 0) {
3141 ijk[0] += LeafNode::DIM;
3142 }
else if (VoxelEdgeAcc::AXIS == 1) {
3143 ijk[1] += LeafNode::DIM;
3144 lhsOffsets = &voxels.maxY();
3145 rhsOffsets = &voxels.minY();
3146 }
else if (VoxelEdgeAcc::AXIS == 2) {
3147 ijk[2] += LeafNode::DIM;
3148 lhsOffsets = &voxels.maxZ();
3149 rhsOffsets = &voxels.minZ();
3152 typename LeafNode::ValueType
value;
3153 const LeafNode* rhsNodePt = acc.probeConstLeaf(ijk);
3156 for (
size_t n = 0, N = lhsOffsets->size();
n <
N; ++
n) {
3157 const Index& pos = (*lhsOffsets)[
n];
3158 bool isActive = lhsNode.isValueOn(pos) || rhsNodePt->isValueOn((*rhsOffsets)[
n]);
3159 if (isActive && (isInsideValue(lhsNode.getValue(pos), iso) !=
3160 isInsideValue(rhsNodePt->getValue((*rhsOffsets)[n]), iso))) {
3161 edgeAcc.set(lhsNode.offsetToGlobalCoord(pos));
3164 }
else if (!acc.probeValue(ijk, value)) {
3165 const bool inside = isInsideValue(value, iso);
3166 for (
size_t n = 0, N = lhsOffsets->size(); n <
N; ++
n) {
3167 const Index& pos = (*lhsOffsets)[
n];
3168 if (lhsNode.isValueOn(pos) && (inside != isInsideValue(lhsNode.getValue(pos), iso))) {
3169 edgeAcc.set(lhsNode.offsetToGlobalCoord(pos));
3179 template<
typename LeafNode,
typename TreeAcc,
typename VoxelEdgeAcc>
3181 evalExtrenalVoxelEdgesInv(VoxelEdgeAcc& edgeAcc, TreeAcc& acc,
const LeafNode& leafnode,
3182 const LeafNodeVoxelOffsets& voxels,
const typename LeafNode::ValueType iso)
3184 Coord ijk = leafnode.origin();
3185 if (VoxelEdgeAcc::AXIS == 0) --ijk[0];
3186 else if (VoxelEdgeAcc::AXIS == 1) --ijk[1];
3187 else if (VoxelEdgeAcc::AXIS == 2) --ijk[2];
3189 typename LeafNode::ValueType
value;
3190 if (!acc.probeConstLeaf(ijk) && !acc.probeValue(ijk, value)) {
3192 const std::vector<Index>* offsets = &voxels.internalNeighborsX();
3193 if (VoxelEdgeAcc::AXIS == 1) offsets = &voxels.internalNeighborsY();
3194 else if (VoxelEdgeAcc::AXIS == 2) offsets = &voxels.internalNeighborsZ();
3196 const bool inside = isInsideValue(value, iso);
3197 for (
size_t n = 0, N = offsets->size(); n <
N; ++
n) {
3199 const Index& pos = (*offsets)[
n];
3200 if (leafnode.isValueOn(pos)
3201 && (inside != isInsideValue(leafnode.getValue(pos), iso)))
3203 ijk = leafnode.offsetToGlobalCoord(pos);
3204 if (VoxelEdgeAcc::AXIS == 0) --ijk[0];
3205 else if (VoxelEdgeAcc::AXIS == 1) --ijk[1];
3206 else if (VoxelEdgeAcc::AXIS == 2) --ijk[2];
3216 template<
typename InputTreeType>
3217 struct IdentifyIntersectingVoxels
3219 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
3220 using InputValueType =
typename InputLeafNodeType::ValueType;
3222 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
3224 IdentifyIntersectingVoxels(
3225 const InputTreeType& inputTree,
3226 const std::vector<const InputLeafNodeType*>& inputLeafNodes,
3227 BoolTreeType& intersectionTree,
3228 InputValueType iso);
3230 IdentifyIntersectingVoxels(IdentifyIntersectingVoxels&,
tbb::split);
3231 void operator()(
const tbb::blocked_range<size_t>&);
3232 void join(
const IdentifyIntersectingVoxels& rhs) {
3233 mIntersectionAccessor.tree().merge(rhs.mIntersectionAccessor.tree());
3237 tree::ValueAccessor<const InputTreeType> mInputAccessor;
3238 InputLeafNodeType
const *
const *
const mInputNodes;
3240 BoolTreeType mIntersectionTree;
3241 tree::ValueAccessor<BoolTreeType> mIntersectionAccessor;
3243 LeafNodeVoxelOffsets mOffsetData;
3244 const LeafNodeVoxelOffsets* mOffsets;
3246 InputValueType mIsovalue;
3250 template<
typename InputTreeType>
3251 IdentifyIntersectingVoxels<InputTreeType>::IdentifyIntersectingVoxels(
3252 const InputTreeType& inputTree,
3253 const std::vector<const InputLeafNodeType*>& inputLeafNodes,
3254 BoolTreeType& intersectionTree,
3256 : mInputAccessor(inputTree)
3257 , mInputNodes(inputLeafNodes.data())
3258 , mIntersectionTree(false)
3259 , mIntersectionAccessor(intersectionTree)
3261 , mOffsets(&mOffsetData)
3264 mOffsetData.constructOffsetList<InputLeafNodeType>();
3268 template<
typename InputTreeType>
3269 IdentifyIntersectingVoxels<InputTreeType>::IdentifyIntersectingVoxels(
3271 : mInputAccessor(rhs.mInputAccessor.tree())
3272 , mInputNodes(rhs.mInputNodes)
3273 , mIntersectionTree(false)
3274 , mIntersectionAccessor(mIntersectionTree)
3276 , mOffsets(rhs.mOffsets)
3277 , mIsovalue(rhs.mIsovalue)
3282 template<
typename InputTreeType>
3284 IdentifyIntersectingVoxels<InputTreeType>::operator()(
const tbb::blocked_range<size_t>& range)
3286 VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 0> xEdgeAcc(mIntersectionAccessor);
3287 VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 1> yEdgeAcc(mIntersectionAccessor);
3288 VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 2> zEdgeAcc(mIntersectionAccessor);
3290 for (
size_t n = range.begin(); n != range.end(); ++
n) {
3292 const InputLeafNodeType& node = *mInputNodes[
n];
3295 evalInternalVoxelEdges(xEdgeAcc, node, *mOffsets, mIsovalue);
3297 evalInternalVoxelEdges(yEdgeAcc, node, *mOffsets, mIsovalue);
3299 evalInternalVoxelEdges(zEdgeAcc, node, *mOffsets, mIsovalue);
3302 evalExtrenalVoxelEdges(xEdgeAcc, mInputAccessor, node, *mOffsets, mIsovalue);
3304 evalExtrenalVoxelEdges(yEdgeAcc, mInputAccessor, node, *mOffsets, mIsovalue);
3306 evalExtrenalVoxelEdges(zEdgeAcc, mInputAccessor, node, *mOffsets, mIsovalue);
3312 evalExtrenalVoxelEdgesInv(xEdgeAcc, mInputAccessor, node, *mOffsets, mIsovalue);
3314 evalExtrenalVoxelEdgesInv(yEdgeAcc, mInputAccessor, node, *mOffsets, mIsovalue);
3316 evalExtrenalVoxelEdgesInv(zEdgeAcc, mInputAccessor, node, *mOffsets, mIsovalue);
3321 template<
typename InputTreeType>
3323 identifySurfaceIntersectingVoxels(
3324 typename InputTreeType::template ValueConverter<bool>::Type& intersectionTree,
3325 const InputTreeType& inputTree,
3326 typename InputTreeType::ValueType isovalue)
3328 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
3330 std::vector<const InputLeafNodeType*> inputLeafNodes;
3331 inputTree.getNodes(inputLeafNodes);
3333 IdentifyIntersectingVoxels<InputTreeType> op(
3334 inputTree, inputLeafNodes, intersectionTree, isovalue);
3336 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, inputLeafNodes.size()), op);
3338 maskActiveTileBorders(inputTree, isovalue, intersectionTree);
3345 template<
typename InputTreeType>
3346 struct MaskIntersectingVoxels
3348 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
3349 using InputValueType =
typename InputLeafNodeType::ValueType;
3351 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
3352 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
3354 MaskIntersectingVoxels(
3355 const InputTreeType& inputTree,
3356 const std::vector<BoolLeafNodeType*>& nodes,
3357 BoolTreeType& intersectionTree,
3358 InputValueType iso);
3360 MaskIntersectingVoxels(MaskIntersectingVoxels&,
tbb::split);
3361 void operator()(
const tbb::blocked_range<size_t>&);
3362 void join(
const MaskIntersectingVoxels& rhs) {
3363 mIntersectionAccessor.tree().merge(rhs.mIntersectionAccessor.tree());
3367 tree::ValueAccessor<const InputTreeType> mInputAccessor;
3368 BoolLeafNodeType
const *
const *
const mNodes;
3370 BoolTreeType mIntersectionTree;
3371 tree::ValueAccessor<BoolTreeType> mIntersectionAccessor;
3373 InputValueType mIsovalue;
3377 template<
typename InputTreeType>
3378 MaskIntersectingVoxels<InputTreeType>::MaskIntersectingVoxels(
3379 const InputTreeType& inputTree,
3380 const std::vector<BoolLeafNodeType*>& nodes,
3381 BoolTreeType& intersectionTree,
3383 : mInputAccessor(inputTree)
3384 , mNodes(nodes.data())
3385 , mIntersectionTree(false)
3386 , mIntersectionAccessor(intersectionTree)
3392 template<
typename InputTreeType>
3393 MaskIntersectingVoxels<InputTreeType>::MaskIntersectingVoxels(
3395 : mInputAccessor(rhs.mInputAccessor.tree())
3396 , mNodes(rhs.mNodes)
3397 , mIntersectionTree(false)
3398 , mIntersectionAccessor(mIntersectionTree)
3399 , mIsovalue(rhs.mIsovalue)
3404 template<
typename InputTreeType>
3406 MaskIntersectingVoxels<InputTreeType>::operator()(
const tbb::blocked_range<size_t>& range)
3408 VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 0> xEdgeAcc(mIntersectionAccessor);
3409 VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 1> yEdgeAcc(mIntersectionAccessor);
3410 VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 2> zEdgeAcc(mIntersectionAccessor);
3413 InputValueType iso(mIsovalue);
3415 for (
size_t n = range.begin(); n != range.end(); ++
n) {
3417 const BoolLeafNodeType& node = *mNodes[
n];
3419 for (
typename BoolLeafNodeType::ValueOnCIter it = node.cbeginValueOn(); it; ++it) {
3421 if (!it.getValue()) {
3423 ijk = it.getCoord();
3425 const bool inside = isInsideValue(mInputAccessor.getValue(ijk), iso);
3427 if (inside != isInsideValue(mInputAccessor.getValue(ijk.offsetBy(1, 0, 0)), iso)) {
3431 if (inside != isInsideValue(mInputAccessor.getValue(ijk.offsetBy(0, 1, 0)), iso)) {
3435 if (inside != isInsideValue(mInputAccessor.getValue(ijk.offsetBy(0, 0, 1)), iso)) {
3444 template<
typename BoolTreeType>
3445 struct MaskBorderVoxels
3447 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
3449 MaskBorderVoxels(
const BoolTreeType& maskTree,
3450 const std::vector<BoolLeafNodeType*>& maskNodes,
3451 BoolTreeType& borderTree)
3452 : mMaskTree(&maskTree)
3453 , mMaskNodes(maskNodes.data())
3454 , mTmpBorderTree(false)
3455 , mBorderTree(&borderTree)
3459 MaskBorderVoxels(MaskBorderVoxels& rhs,
tbb::split)
3460 : mMaskTree(rhs.mMaskTree)
3461 , mMaskNodes(rhs.mMaskNodes)
3462 , mTmpBorderTree(false)
3463 , mBorderTree(&mTmpBorderTree)
3467 void join(MaskBorderVoxels& rhs) { mBorderTree->merge(*rhs.mBorderTree); }
3469 void operator()(
const tbb::blocked_range<size_t>& range)
3471 tree::ValueAccessor<const BoolTreeType> maskAcc(*mMaskTree);
3472 tree::ValueAccessor<BoolTreeType> borderAcc(*mBorderTree);
3475 for (
size_t n = range.begin(); n != range.end(); ++
n) {
3477 const BoolLeafNodeType& node = *mMaskNodes[
n];
3479 for (
typename BoolLeafNodeType::ValueOnCIter it = node.cbeginValueOn(); it; ++it) {
3481 ijk = it.getCoord();
3483 const bool lhs = it.getValue();
3486 bool isEdgeVoxel =
false;
3489 isEdgeVoxel = (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3492 isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3495 isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3498 isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3502 isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3505 isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3508 isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3512 borderAcc.setValue(ijk,
true);
3519 BoolTreeType
const *
const mMaskTree;
3520 BoolLeafNodeType
const *
const *
const mMaskNodes;
3522 BoolTreeType mTmpBorderTree;
3523 BoolTreeType *
const mBorderTree;
3527 template<
typename BoolTreeType>
3528 struct SyncMaskValues
3530 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
3532 SyncMaskValues(
const std::vector<BoolLeafNodeType*>& nodes,
const BoolTreeType& mask)
3533 : mNodes(nodes.data())
3538 void operator()(
const tbb::blocked_range<size_t>& range)
const
3540 using ValueOnIter =
typename BoolLeafNodeType::ValueOnIter;
3542 tree::ValueAccessor<const BoolTreeType> maskTreeAcc(*mMaskTree);
3544 for (
size_t n = range.begin(), N = range.end(); n !=
N; ++
n) {
3546 BoolLeafNodeType& node = *mNodes[
n];
3548 const BoolLeafNodeType * maskNode = maskTreeAcc.probeConstLeaf(node.origin());
3551 for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3552 const Index pos = it.pos();
3553 if (maskNode->getValue(pos)) {
3554 node.setValueOnly(pos,
true);
3562 BoolLeafNodeType *
const *
const mNodes;
3563 BoolTreeType
const *
const mMaskTree;
3570 template<
typename BoolTreeType>
3573 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
3575 MaskSurface(
const std::vector<BoolLeafNodeType*>& nodes,
const BoolTreeType& mask,
3576 const math::Transform& inputTransform,
const math::Transform& maskTransform,
bool invert)
3577 : mNodes(nodes.data())
3579 , mInputTransform(inputTransform)
3580 , mMaskTransform(maskTransform)
3581 , mInvertMask(invert)
3585 void operator()(
const tbb::blocked_range<size_t>& range)
const
3587 using ValueOnIter =
typename BoolLeafNodeType::ValueOnIter;
3589 tree::ValueAccessor<const BoolTreeType> maskTreeAcc(*mMaskTree);
3591 const bool matchingTransforms = mInputTransform == mMaskTransform;
3592 const bool maskState = mInvertMask;
3594 for (
size_t n = range.begin(), N = range.end(); n !=
N; ++
n) {
3596 BoolLeafNodeType& node = *mNodes[
n];
3598 if (matchingTransforms) {
3600 const BoolLeafNodeType * maskNode = maskTreeAcc.probeConstLeaf(node.origin());
3604 for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3605 const Index pos = it.pos();
3606 if (maskNode->isValueOn(pos) == maskState) {
3607 node.setValueOnly(pos,
true);
3613 if (maskTreeAcc.isValueOn(node.origin()) == maskState) {
3614 for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3615 node.setValueOnly(it.pos(),
true);
3625 for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3627 ijk = mMaskTransform.worldToIndexCellCentered(
3628 mInputTransform.indexToWorld(it.getCoord()));
3630 if (maskTreeAcc.isValueOn(ijk) == maskState) {
3631 node.setValueOnly(it.pos(),
true);
3640 BoolLeafNodeType *
const *
const mNodes;
3641 BoolTreeType
const *
const mMaskTree;
3642 math::Transform
const mInputTransform;
3643 math::Transform
const mMaskTransform;
3644 bool const mInvertMask;
3648 template<
typename InputGr
idType>
3651 typename InputGridType::TreeType::template ValueConverter<bool>::Type& intersectionTree,
3652 typename InputGridType::TreeType::template ValueConverter<bool>::Type& borderTree,
3653 const InputGridType& inputGrid,
3656 typename InputGridType::ValueType isovalue)
3658 using InputTreeType =
typename InputGridType::TreeType;
3659 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
3660 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
3661 using BoolGridType = Grid<BoolTreeType>;
3663 if (maskGrid && maskGrid->type() == BoolGridType::gridType()) {
3665 const math::Transform&
transform = inputGrid.transform();
3666 const InputTreeType& inputTree = inputGrid.tree();
3668 const BoolGridType * surfaceMask =
static_cast<const BoolGridType*
>(maskGrid.get());
3670 const BoolTreeType& maskTree = surfaceMask->tree();
3671 const math::Transform& maskTransform = surfaceMask->transform();
3676 std::vector<BoolLeafNodeType*> intersectionLeafNodes;
3677 intersectionTree.getNodes(intersectionLeafNodes);
3680 MaskSurface<BoolTreeType>(
3681 intersectionLeafNodes, maskTree, transform, maskTransform, invertMask));
3686 MaskBorderVoxels<BoolTreeType> borderOp(
3687 intersectionTree, intersectionLeafNodes, borderTree);
3688 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, intersectionLeafNodes.size()), borderOp);
3693 BoolTreeType tmpIntersectionTree(
false);
3695 MaskIntersectingVoxels<InputTreeType> op(
3696 inputTree, intersectionLeafNodes, tmpIntersectionTree, isovalue);
3698 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, intersectionLeafNodes.size()), op);
3700 std::vector<BoolLeafNodeType*> tmpIntersectionLeafNodes;
3701 tmpIntersectionTree.getNodes(tmpIntersectionLeafNodes);
3703 tbb::parallel_for(tbb::blocked_range<size_t>(0, tmpIntersectionLeafNodes.size()),
3704 SyncMaskValues<BoolTreeType>(tmpIntersectionLeafNodes, intersectionTree));
3706 intersectionTree.clear();
3707 intersectionTree.merge(tmpIntersectionTree);
3715 template<
typename InputTreeType>
3716 struct ComputeAuxiliaryData
3718 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
3719 using InputValueType =
typename InputLeafNodeType::ValueType;
3721 using BoolLeafNodeType = tree::LeafNode<bool, InputLeafNodeType::LOG2DIM>;
3727 ComputeAuxiliaryData(
const InputTreeType& inputTree,
3728 const std::vector<const BoolLeafNodeType*>& intersectionLeafNodes,
3729 Int16TreeType& signFlagsTree,
3730 Index32TreeType& pointIndexTree,
3731 InputValueType iso);
3733 ComputeAuxiliaryData(ComputeAuxiliaryData&,
tbb::split);
3734 void operator()(
const tbb::blocked_range<size_t>&);
3735 void join(
const ComputeAuxiliaryData& rhs) {
3736 mSignFlagsAccessor.tree().merge(rhs.mSignFlagsAccessor.tree());
3737 mPointIndexAccessor.tree().merge(rhs.mPointIndexAccessor.tree());
3741 tree::ValueAccessor<const InputTreeType> mInputAccessor;
3742 BoolLeafNodeType
const *
const *
const mIntersectionNodes;
3744 Int16TreeType mSignFlagsTree;
3745 tree::ValueAccessor<Int16TreeType> mSignFlagsAccessor;
3746 Index32TreeType mPointIndexTree;
3747 tree::ValueAccessor<Index32TreeType> mPointIndexAccessor;
3749 const InputValueType mIsovalue;
3753 template<
typename InputTreeType>
3754 ComputeAuxiliaryData<InputTreeType>::ComputeAuxiliaryData(
3755 const InputTreeType& inputTree,
3756 const std::vector<const BoolLeafNodeType*>& intersectionLeafNodes,
3757 Int16TreeType& signFlagsTree,
3758 Index32TreeType& pointIndexTree,
3760 : mInputAccessor(inputTree)
3761 , mIntersectionNodes(intersectionLeafNodes.data())
3763 , mSignFlagsAccessor(signFlagsTree)
3764 , mPointIndexTree(std::numeric_limits<
Index32>::
max())
3765 , mPointIndexAccessor(pointIndexTree)
3772 template<
typename InputTreeType>
3773 ComputeAuxiliaryData<InputTreeType>::ComputeAuxiliaryData(ComputeAuxiliaryData& rhs,
tbb::split)
3774 : mInputAccessor(rhs.mInputAccessor.tree())
3775 , mIntersectionNodes(rhs.mIntersectionNodes)
3777 , mSignFlagsAccessor(mSignFlagsTree)
3778 , mPointIndexTree(std::numeric_limits<
Index32>::
max())
3779 , mPointIndexAccessor(mPointIndexTree)
3780 , mIsovalue(rhs.mIsovalue)
3785 template<
typename InputTreeType>
3787 ComputeAuxiliaryData<InputTreeType>::operator()(
const tbb::blocked_range<size_t>& range)
3789 using Int16LeafNodeType =
typename Int16TreeType::LeafNodeType;
3792 math::Tuple<8, InputValueType> cellVertexValues;
3793 typename std::unique_ptr<Int16LeafNodeType> signsNodePt(
new Int16LeafNodeType(ijk, 0));
3795 for (
size_t n = range.begin(), N = range.end(); n !=
N; ++
n) {
3797 const BoolLeafNodeType& maskNode = *mIntersectionNodes[
n];
3798 const Coord& origin = maskNode.origin();
3800 const InputLeafNodeType *leafPt = mInputAccessor.probeConstLeaf(origin);
3802 if (!signsNodePt.get()) signsNodePt.reset(
new Int16LeafNodeType(origin, 0));
3803 else signsNodePt->setOrigin(origin);
3805 bool updatedNode =
false;
3807 for (
typename BoolLeafNodeType::ValueOnCIter it = maskNode.cbeginValueOn(); it; ++it) {
3809 const Index pos = it.pos();
3810 ijk = BoolLeafNodeType::offsetToLocalCoord(pos);
3813 ijk[0] <
int(BoolLeafNodeType::DIM - 1) &&
3814 ijk[1] <
int(BoolLeafNodeType::DIM - 1) &&
3815 ijk[2] <
int(BoolLeafNodeType::DIM - 1) ) {
3816 getCellVertexValues(*leafPt, pos, cellVertexValues);
3818 getCellVertexValues(mInputAccessor, origin + ijk, cellVertexValues);
3821 uint8_t signFlags = computeSignFlags(cellVertexValues, mIsovalue);
3823 if (signFlags != 0 && signFlags != 0xFF) {
3825 const bool inside = signFlags & 0x1;
3827 int edgeFlags = inside ? INSIDE : 0;
3829 if (!it.getValue()) {
3830 edgeFlags |= inside != ((signFlags & 0x02) != 0) ? XEDGE : 0;
3831 edgeFlags |= inside != ((signFlags & 0x10) != 0) ? YEDGE : 0;
3832 edgeFlags |= inside != ((signFlags & 0x08) != 0) ? ZEDGE : 0;
3835 const uint8_t ambiguousCellFlags = sAmbiguousFace[signFlags];
3836 if (ambiguousCellFlags != 0) {
3837 correctCellSigns(signFlags, ambiguousCellFlags, mInputAccessor,
3838 origin + ijk, mIsovalue);
3841 edgeFlags |=
int(signFlags);
3843 signsNodePt->setValueOn(pos,
Int16(edgeFlags));
3849 typename Index32TreeType::LeafNodeType* idxNode = mPointIndexAccessor.touchLeaf(origin);
3850 idxNode->topologyUnion(*signsNodePt);
3853 for (
auto it = idxNode->beginValueOn(); it; ++it) {
3854 idxNode->setValueOnly(it.pos(), 0);
3857 mSignFlagsAccessor.addLeaf(signsNodePt.release());
3863 template<
typename InputTreeType>
3865 computeAuxiliaryData(
3868 const typename InputTreeType::template ValueConverter<bool>::Type& intersectionTree,
3869 const InputTreeType& inputTree,
3870 typename InputTreeType::ValueType isovalue)
3872 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
3873 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
3875 std::vector<const BoolLeafNodeType*> intersectionLeafNodes;
3876 intersectionTree.getNodes(intersectionLeafNodes);
3878 ComputeAuxiliaryData<InputTreeType> op(
3879 inputTree, intersectionLeafNodes, signFlagsTree, pointIndexTree, isovalue);
3881 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, intersectionLeafNodes.size()), op);
3888 template<Index32 LeafNodeLog2Dim>
3889 struct LeafNodePointCount
3891 using Int16LeafNodeType = tree::LeafNode<Int16, LeafNodeLog2Dim>;
3893 LeafNodePointCount(
const std::vector<Int16LeafNodeType*>& leafNodes,
3894 std::unique_ptr<
Index32[]>& leafNodeCount)
3895 : mLeafNodes(leafNodes.data())
3896 , mData(leafNodeCount.get())
3900 void operator()(
const tbb::blocked_range<size_t>& range)
const {
3902 for (
size_t n = range.begin(), N = range.end(); n !=
N; ++
n) {
3906 Int16 const * p = mLeafNodes[
n]->buffer().data();
3910 count +=
Index32(sEdgeGroupTable[(SIGNS & *p)][0]);
3919 Int16LeafNodeType *
const *
const mLeafNodes;
3924 template<
typename Po
intIndexLeafNode>
3925 struct AdaptiveLeafNodePointCount
3927 using Int16LeafNodeType = tree::LeafNode<Int16, PointIndexLeafNode::LOG2DIM>;
3929 AdaptiveLeafNodePointCount(
const std::vector<PointIndexLeafNode*>& pointIndexNodes,
3930 const std::vector<Int16LeafNodeType*>& signDataNodes,
3931 std::unique_ptr<
Index32[]>& leafNodeCount)
3932 : mPointIndexNodes(pointIndexNodes.data())
3933 , mSignDataNodes(signDataNodes.data())
3934 , mData(leafNodeCount.get())
3938 void operator()(
const tbb::blocked_range<size_t>& range)
const
3940 using IndexType =
typename PointIndexLeafNode::ValueType;
3942 for (
size_t n = range.begin(), N = range.end(); n !=
N; ++
n) {
3944 const PointIndexLeafNode& node = *mPointIndexNodes[
n];
3945 const Int16LeafNodeType& signNode = *mSignDataNodes[
n];
3949 std::set<IndexType> uniqueRegions;
3951 for (
typename PointIndexLeafNode::ValueOnCIter it = node.cbeginValueOn(); it; ++it) {
3953 IndexType
id = it.getValue();
3956 count += size_t(sEdgeGroupTable[(SIGNS & signNode.getValue(it.pos()))][0]);
3958 uniqueRegions.insert(
id);
3962 mData[
n] =
Index32(count + uniqueRegions.size());
3967 PointIndexLeafNode
const *
const *
const mPointIndexNodes;
3968 Int16LeafNodeType
const *
const *
const mSignDataNodes;
3973 template<
typename Po
intIndexLeafNode>
3976 using Int16LeafNodeType = tree::LeafNode<Int16, PointIndexLeafNode::LOG2DIM>;
3978 MapPoints(
const std::vector<PointIndexLeafNode*>& pointIndexNodes,
3979 const std::vector<Int16LeafNodeType*>& signDataNodes,
3980 std::unique_ptr<
Index32[]>& leafNodeCount)
3981 : mPointIndexNodes(pointIndexNodes.data())
3982 , mSignDataNodes(signDataNodes.data())
3983 , mData(leafNodeCount.get())
3987 void operator()(
const tbb::blocked_range<size_t>& range)
const {
3989 for (
size_t n = range.begin(), N = range.end(); n !=
N; ++
n) {
3991 const Int16LeafNodeType& signNode = *mSignDataNodes[
n];
3992 PointIndexLeafNode& indexNode = *mPointIndexNodes[
n];
3996 for (
auto it = indexNode.beginValueOn(); it; ++it) {
3997 const Index pos = it.pos();
3998 indexNode.setValueOnly(pos, pointOffset);
3999 const int signs = SIGNS &
int(signNode.getValue(pos));
4000 pointOffset +=
Index32(sEdgeGroupTable[signs][0]);
4006 PointIndexLeafNode *
const *
const mPointIndexNodes;
4007 Int16LeafNodeType
const *
const *
const mSignDataNodes;
4014 template<
typename TreeType,
typename PrimBuilder>
4015 struct ComputePolygons
4018 using Int16LeafNodeType =
typename Int16TreeType::LeafNodeType;
4021 using Index32LeafNodeType =
typename Index32TreeType::LeafNodeType;
4025 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
4026 const Int16TreeType& signFlagsTree,
4027 const Index32TreeType& idxTree,
4029 bool invertSurfaceOrientation);
4031 void setRefSignTree(
const Int16TreeType *
r) { mRefSignFlagsTree =
r; }
4033 void operator()(
const tbb::blocked_range<size_t>&)
const;
4036 Int16LeafNodeType *
const *
const mSignFlagsLeafNodes;
4037 Int16TreeType
const *
const mSignFlagsTree;
4038 Int16TreeType
const * mRefSignFlagsTree;
4039 Index32TreeType
const *
const mIndexTree;
4041 bool const mInvertSurfaceOrientation;
4045 template<
typename TreeType,
typename PrimBuilder>
4046 ComputePolygons<TreeType, PrimBuilder>::ComputePolygons(
4047 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
4048 const Int16TreeType& signFlagsTree,
4049 const Index32TreeType& idxTree,
4051 bool invertSurfaceOrientation)
4052 : mSignFlagsLeafNodes(signFlagsLeafNodes.data())
4053 , mSignFlagsTree(&signFlagsTree)
4054 , mRefSignFlagsTree(nullptr)
4055 , mIndexTree(&idxTree)
4056 , mPolygonPoolList(&polygons)
4057 , mInvertSurfaceOrientation(invertSurfaceOrientation)
4061 template<
typename InputTreeType,
typename PrimBuilder>
4063 ComputePolygons<InputTreeType, PrimBuilder>::operator()(
const tbb::blocked_range<size_t>& range)
const
4065 using Int16ValueAccessor = tree::ValueAccessor<const Int16TreeType>;
4066 Int16ValueAccessor signAcc(*mSignFlagsTree);
4068 tree::ValueAccessor<const Index32TreeType> idxAcc(*mIndexTree);
4070 const bool invertSurfaceOrientation = mInvertSurfaceOrientation;
4077 std::unique_ptr<Int16ValueAccessor> refSignAcc;
4078 if (mRefSignFlagsTree) refSignAcc.reset(
new Int16ValueAccessor(*mRefSignFlagsTree));
4080 for (
size_t n = range.begin(); n != range.end(); ++
n) {
4082 const Int16LeafNodeType& node = *mSignFlagsLeafNodes[
n];
4083 origin = node.origin();
4087 typename Int16LeafNodeType::ValueOnCIter iter = node.cbeginValueOn();
4088 for (; iter; ++iter) {
4089 if (iter.getValue() & XEDGE) ++edgeCount;
4090 if (iter.getValue() & YEDGE) ++edgeCount;
4091 if (iter.getValue() & ZEDGE) ++edgeCount;
4094 if(edgeCount == 0)
continue;
4096 mesher.init(edgeCount, (*mPolygonPoolList)[n]);
4098 const Int16LeafNodeType *signleafPt = signAcc.probeConstLeaf(origin);
4099 const Index32LeafNodeType *idxLeafPt = idxAcc.probeConstLeaf(origin);
4101 if (!signleafPt || !idxLeafPt)
continue;
4104 const Int16LeafNodeType *refSignLeafPt =
nullptr;
4105 if (refSignAcc) refSignLeafPt = refSignAcc->probeConstLeaf(origin);
4109 for (iter = node.cbeginValueOn(); iter; ++iter) {
4110 ijk = iter.getCoord();
4112 Int16 flags = iter.getValue();
4114 if (!(flags & 0xE00))
continue;
4117 if (refSignLeafPt) {
4118 refFlags = refSignLeafPt->getValue(iter.pos());
4125 const uint8_t cell = uint8_t(SIGNS & flags);
4127 if (sEdgeGroupTable[cell][0] > 1) {
4128 offsets[0] = (sEdgeGroupTable[cell][1] - 1);
4129 offsets[1] = (sEdgeGroupTable[cell][9] - 1);
4130 offsets[2] = (sEdgeGroupTable[cell][4] - 1);
4133 if (ijk[0] > origin[0] && ijk[1] > origin[1] && ijk[2] > origin[2]) {
4134 constructPolygons(invertSurfaceOrientation,
4135 flags, refFlags, offsets, ijk, *signleafPt, *idxLeafPt, mesher);
4137 constructPolygons(invertSurfaceOrientation,
4138 flags, refFlags, offsets, ijk, signAcc, idxAcc, mesher);
4151 template<
typename T>
4154 CopyArray(T * outputArray,
const T * inputArray,
size_t outputOffset = 0)
4155 : mOutputArray(outputArray), mInputArray(inputArray), mOutputOffset(outputOffset)
4159 void operator()(
const tbb::blocked_range<size_t>& inputArrayRange)
const
4161 const size_t offset = mOutputOffset;
4162 for (
size_t n = inputArrayRange.begin(), N = inputArrayRange.end(); n <
N; ++
n) {
4163 mOutputArray[offset +
n] = mInputArray[
n];
4168 T *
const mOutputArray;
4169 T const *
const mInputArray;
4170 size_t const mOutputOffset;
4174 struct FlagAndCountQuadsToSubdivide
4177 const std::vector<uint8_t>& pointFlags,
4179 std::unique_ptr<
unsigned[]>& numQuadsToDivide)
4180 : mPolygonPoolList(&polygons)
4181 , mPointFlags(pointFlags.data())
4182 , mPoints(points.get())
4183 , mNumQuadsToDivide(numQuadsToDivide.get())
4187 void operator()(
const tbb::blocked_range<size_t>& range)
const
4189 for (
size_t n = range.begin(), N = range.end(); n <
N; ++
n) {
4191 PolygonPool& polygons = (*mPolygonPoolList)[
n];
4196 for (
size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4198 char& flags = polygons.quadFlags(i);
4202 Vec4I& quad = polygons.quad(i);
4204 const bool edgePoly = mPointFlags[quad[0]] || mPointFlags[quad[1]]
4205 || mPointFlags[quad[2]] || mPointFlags[quad[3]];
4207 if (!edgePoly)
continue;
4209 const Vec3s& p0 = mPoints[quad[0]];
4210 const Vec3s& p1 = mPoints[quad[1]];
4211 const Vec3s& p2 = mPoints[quad[2]];
4212 const Vec3s& p3 = mPoints[quad[3]];
4214 if (!isPlanarQuad(p0, p1, p2, p3, 1e-6
f)) {
4221 mNumQuadsToDivide[
n] =
count;
4227 uint8_t
const *
const mPointFlags;
4228 Vec3s const *
const mPoints;
4229 unsigned *
const mNumQuadsToDivide;
4233 struct SubdivideQuads
4239 std::unique_ptr<
unsigned[]>& numQuadsToDivide,
4240 std::unique_ptr<
unsigned[]>& centroidOffsets)
4241 : mPolygonPoolList(&polygons)
4242 , mPoints(points.get())
4243 , mCentroids(centroids.get())
4244 , mNumQuadsToDivide(numQuadsToDivide.get())
4245 , mCentroidOffsets(centroidOffsets.get())
4246 , mPointCount(pointCount)
4250 void operator()(
const tbb::blocked_range<size_t>& range)
const
4252 for (
size_t n = range.begin(), N = range.end(); n <
N; ++
n) {
4254 PolygonPool& polygons = (*mPolygonPoolList)[
n];
4256 const size_t nonplanarCount = size_t(mNumQuadsToDivide[n]);
4258 if (nonplanarCount > 0) {
4260 PolygonPool tmpPolygons;
4261 tmpPolygons.resetQuads(polygons.numQuads() - nonplanarCount);
4262 tmpPolygons.resetTriangles(polygons.numTriangles() + size_t(4) * nonplanarCount);
4264 size_t offset = mCentroidOffsets[
n];
4266 size_t triangleIdx = 0;
4268 for (
size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4270 const char quadFlags = polygons.quadFlags(i);
4273 unsigned newPointIdx = unsigned(offset + mPointCount);
4277 mCentroids[
offset] = (mPoints[quad[0]] + mPoints[quad[1]] +
4278 mPoints[quad[2]] + mPoints[quad[3]]) * 0.25
f;
4283 Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
4285 triangle[0] = quad[0];
4286 triangle[1] = newPointIdx;
4287 triangle[2] = quad[3];
4289 tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
4295 Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
4297 triangle[0] = quad[0];
4298 triangle[1] = quad[1];
4299 triangle[2] = newPointIdx;
4301 tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
4307 Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
4309 triangle[0] = quad[1];
4310 triangle[1] = quad[2];
4311 triangle[2] = newPointIdx;
4313 tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
4320 Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
4322 triangle[0] = quad[2];
4323 triangle[1] = quad[3];
4324 triangle[2] = newPointIdx;
4326 tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
4335 for (
size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
4336 tmpPolygons.triangle(triangleIdx) = polygons.triangle(i);
4337 tmpPolygons.triangleFlags(triangleIdx) = polygons.triangleFlags(i);
4342 for (
size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4346 tmpPolygons.quad(quadIdx) = quad;
4347 tmpPolygons.quadFlags(quadIdx) = polygons.quadFlags(i);
4352 polygons.copy(tmpPolygons);
4359 Vec3s const *
const mPoints;
4360 Vec3s *
const mCentroids;
4361 unsigned *
const mNumQuadsToDivide;
4362 unsigned *
const mCentroidOffsets;
4363 size_t const mPointCount;
4368 subdivideNonplanarSeamLineQuads(
4370 size_t polygonPoolListSize,
4372 size_t& pointListSize,
4373 std::vector<uint8_t>& pointFlags)
4375 const tbb::blocked_range<size_t> polygonPoolListRange(0, polygonPoolListSize);
4377 std::unique_ptr<unsigned[]> numQuadsToDivide(
new unsigned[polygonPoolListSize]);
4380 FlagAndCountQuadsToSubdivide(polygonPoolList, pointFlags, pointList, numQuadsToDivide));
4382 std::unique_ptr<unsigned[]> centroidOffsets(
new unsigned[polygonPoolListSize]);
4384 size_t centroidCount = 0;
4388 for (
size_t n = 0, N = polygonPoolListSize; n <
N; ++
n) {
4389 centroidOffsets[
n] = sum;
4390 sum += numQuadsToDivide[
n];
4392 centroidCount = size_t(sum);
4395 std::unique_ptr<Vec3s[]> centroidList(
new Vec3s[centroidCount]);
4398 SubdivideQuads(polygonPoolList, pointList, pointListSize,
4399 centroidList, numQuadsToDivide, centroidOffsets));
4401 if (centroidCount > 0) {
4403 const size_t newPointListSize = centroidCount + pointListSize;
4405 std::unique_ptr<openvdb::Vec3s[]> newPointList(
new openvdb::Vec3s[newPointListSize]);
4408 CopyArray<Vec3s>(newPointList.get(), pointList.get()));
4410 tbb::parallel_for(tbb::blocked_range<size_t>(0, newPointListSize - pointListSize),
4411 CopyArray<Vec3s>(newPointList.get(), centroidList.get(), pointListSize));
4413 pointListSize = newPointListSize;
4414 pointList.swap(newPointList);
4415 pointFlags.resize(pointListSize, 0);
4420 struct ReviseSeamLineFlags
4423 const std::vector<uint8_t>& pointFlags)
4424 : mPolygonPoolList(&polygons)
4425 , mPointFlags(pointFlags.data())
4429 void operator()(
const tbb::blocked_range<size_t>& range)
const
4431 for (
size_t n = range.begin(), N = range.end(); n <
N; ++
n) {
4433 PolygonPool& polygons = (*mPolygonPoolList)[
n];
4435 for (
size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4437 char& flags = polygons.quadFlags(i);
4439 if (flags & POLYFLAG_FRACTURE_SEAM) {
4443 const bool hasSeamLinePoint =
4444 mPointFlags[verts[0]] || mPointFlags[verts[1]] ||
4445 mPointFlags[verts[2]] || mPointFlags[verts[3]];
4447 if (!hasSeamLinePoint) {
4448 flags &= ~POLYFLAG_FRACTURE_SEAM;
4453 for (
size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
4455 char& flags = polygons.triangleFlags(i);
4457 if (flags & POLYFLAG_FRACTURE_SEAM) {
4461 const bool hasSeamLinePoint =
4462 mPointFlags[verts[0]] || mPointFlags[verts[1]] || mPointFlags[verts[2]];
4464 if (!hasSeamLinePoint) {
4465 flags &= ~POLYFLAG_FRACTURE_SEAM;
4476 uint8_t
const *
const mPointFlags;
4481 reviseSeamLineFlags(
PolygonPoolList& polygonPoolList,
size_t polygonPoolListSize,
4482 std::vector<uint8_t>& pointFlags)
4485 ReviseSeamLineFlags(polygonPoolList, pointFlags));
4492 template<
typename InputTreeType>
4493 struct MaskDisorientedTrianglePoints
4495 MaskDisorientedTrianglePoints(
const InputTreeType& inputTree,
const PolygonPoolList& polygons,
4496 const PointList& pointList, std::unique_ptr<uint8_t[]>& pointMask,
4497 const math::Transform& transform,
bool invertSurfaceOrientation)
4498 : mInputTree(&inputTree)
4499 , mPolygonPoolList(&polygons)
4500 , mPointList(&pointList)
4501 , mPointMask(pointMask.get())
4502 , mTransform(transform)
4503 , mInvertSurfaceOrientation(invertSurfaceOrientation)
4507 void operator()(
const tbb::blocked_range<size_t>& range)
const
4509 using ValueType =
typename InputTreeType::LeafNodeType::ValueType;
4511 tree::ValueAccessor<const InputTreeType> inputAcc(*mInputTree);
4512 Vec3s centroid, normal;
4515 const bool invertGradientDir = mInvertSurfaceOrientation || isBoolValue<ValueType>();
4517 for (
size_t n = range.begin(), N = range.end(); n <
N; ++
n) {
4519 const PolygonPool& polygons = (*mPolygonPoolList)[
n];
4521 for (
size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
4523 const Vec3I& verts = polygons.triangle(i);
4525 const Vec3s& v0 = (*mPointList)[verts[0]];
4526 const Vec3s& v1 = (*mPointList)[verts[1]];
4527 const Vec3s& v2 = (*mPointList)[verts[2]];
4529 normal = (v2 -
v0).
cross((v1 - v0));
4532 centroid = (v0 + v1 +
v2) * (1.0
f / 3.0
f);
4533 ijk = mTransform.worldToIndexCellCentered(centroid);
4538 if (invertGradientDir) {
4543 if (dir.dot(normal) < -0.5f) {
4548 mPointMask[verts[0]] = 1;
4549 mPointMask[verts[1]] = 1;
4550 mPointMask[verts[2]] = 1;
4559 InputTreeType
const *
const mInputTree;
4562 uint8_t *
const mPointMask;
4563 math::Transform
const mTransform;
4564 bool const mInvertSurfaceOrientation;
4568 template<
typename InputTree>
4570 relaxDisorientedTriangles(
4571 bool invertSurfaceOrientation,
4572 const InputTree& inputTree,
4573 const math::Transform& transform,
4575 size_t polygonPoolListSize,
4577 const size_t pointListSize)
4579 const tbb::blocked_range<size_t> polygonPoolListRange(0, polygonPoolListSize);
4581 std::unique_ptr<uint8_t[]> pointMask(
new uint8_t[pointListSize]);
4582 fillArray(pointMask.get(), uint8_t(0), pointListSize);
4585 MaskDisorientedTrianglePoints<InputTree>(
4586 inputTree, polygonPoolList, pointList, pointMask, transform, invertSurfaceOrientation));
4588 std::unique_ptr<uint8_t[]> pointUpdates(
new uint8_t[pointListSize]);
4589 fillArray(pointUpdates.get(), uint8_t(0), pointListSize);
4591 std::unique_ptr<Vec3s[]> newPoints(
new Vec3s[pointListSize]);
4592 fillArray(newPoints.get(),
Vec3s(0.0
f, 0.0
f, 0.0
f), pointListSize);
4594 for (
size_t n = 0, N = polygonPoolListSize; n <
N; ++
n) {
4596 PolygonPool& polygons = polygonPoolList[
n];
4598 for (
size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4601 for (
int v = 0; v < 4; ++
v) {
4603 const unsigned pointIdx = verts[
v];
4605 if (pointMask[pointIdx] == 1) {
4607 newPoints[pointIdx] +=
4608 pointList[verts[0]] + pointList[verts[1]] +
4609 pointList[verts[2]] + pointList[verts[3]];
4611 pointUpdates[pointIdx] = uint8_t(pointUpdates[pointIdx] + 4);
4616 for (
size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
4619 for (
int v = 0; v < 3; ++
v) {
4621 const unsigned pointIdx = verts[
v];
4623 if (pointMask[pointIdx] == 1) {
4624 newPoints[pointIdx] +=
4625 pointList[verts[0]] + pointList[verts[1]] + pointList[verts[2]];
4627 pointUpdates[pointIdx] = uint8_t(pointUpdates[pointIdx] + 3);
4633 for (
size_t n = 0, N = pointListSize; n <
N; ++
n) {
4634 if (pointUpdates[n] > 0) {
4635 const double weight = 1.0 / double(pointUpdates[n]);
4636 pointList[
n] = newPoints[
n] * float(weight);
4650 PolygonPool::PolygonPool()
4654 , mTriangles(nullptr)
4655 , mQuadFlags(nullptr)
4656 , mTriangleFlags(nullptr)
4663 : mNumQuads(numQuads)
4664 , mNumTriangles(numTriangles)
4665 , mQuads(new openvdb::
Vec4I[mNumQuads])
4666 , mTriangles(new openvdb::
Vec3I[mNumTriangles])
4667 , mQuadFlags(new char[mNumQuads])
4668 , mTriangleFlags(new char[mNumTriangles])
4679 for (
size_t i = 0; i < mNumQuads; ++i) {
4680 mQuads[i] = rhs.mQuads[i];
4681 mQuadFlags[i] = rhs.mQuadFlags[i];
4684 for (
size_t i = 0; i < mNumTriangles; ++i) {
4685 mTriangles[i] = rhs.mTriangles[i];
4686 mTriangleFlags[i] = rhs.mTriangleFlags[i];
4696 mQuadFlags.reset(
new char[mNumQuads]);
4704 mQuads.reset(
nullptr);
4705 mQuadFlags.reset(
nullptr);
4712 mNumTriangles =
size;
4714 mTriangleFlags.reset(
new char[mNumTriangles]);
4722 mTriangles.reset(
nullptr);
4723 mTriangleFlags.reset(
nullptr);
4730 if (!(n < mNumQuads))
return false;
4735 mQuads.reset(
nullptr);
4739 std::unique_ptr<char[]>
flags(
new char[n]);
4741 for (
size_t i = 0; i <
n; ++i) {
4742 quads[i] = mQuads[i];
4743 flags[i] = mQuadFlags[i];
4747 mQuadFlags.swap(flags);
4759 if (!(n < mNumTriangles))
return false;
4764 mTriangles.reset(
nullptr);
4767 std::unique_ptr<openvdb::Vec3I[]> triangles(
new openvdb::Vec3I[n]);
4768 std::unique_ptr<char[]>
flags(
new char[n]);
4770 for (
size_t i = 0; i <
n; ++i) {
4771 triangles[i] = mTriangles[i];
4772 flags[i] = mTriangleFlags[i];
4775 mTriangles.swap(triangles);
4776 mTriangleFlags.swap(flags);
4793 , mSeamPointListSize(0)
4794 , mPolygonPoolListSize(0)
4795 , mIsovalue(isovalue)
4796 , mPrimAdaptivity(adaptivity)
4797 , mSecAdaptivity(0.0)
4799 , mSurfaceMaskGrid(
GridBase::ConstPtr())
4800 , mAdaptivityGrid(
GridBase::ConstPtr())
4801 , mAdaptivityMaskTree(
TreeBase::ConstPtr())
4804 , mInvertSurfaceMask(false)
4805 , mRelaxDisorientedTriangles(relaxDisorientedTriangles)
4806 , mQuantizedSeamPoints(nullptr)
4816 mSecAdaptivity = secAdaptivity;
4821 mSeamPointListSize = 0;
4822 mQuantizedSeamPoints.reset(
nullptr);
4829 mSurfaceMaskGrid =
mask;
4830 mInvertSurfaceMask = invertMask;
4837 mAdaptivityGrid = grid;
4844 mAdaptivityMaskTree = tree;
4848 template<
typename InputGr
idType>
4854 using InputTreeType =
typename InputGridType::TreeType;
4855 using InputLeafNodeType =
typename InputTreeType::LeafNodeType;
4856 using InputValueType =
typename InputLeafNodeType::ValueType;
4862 using BoolTreeType =
typename InputTreeType::template ValueConverter<bool>::Type;
4864 using Int16LeafNodeType =
typename Int16TreeType::LeafNodeType;
4866 using Index32LeafNodeType =
typename Index32TreeType::LeafNodeType;
4871 mPolygonPoolListSize = 0;
4873 mPointFlags.clear();
4878 const InputValueType isovalue = InputValueType(mIsovalue);
4879 const float adaptivityThreshold = float(mPrimAdaptivity);
4880 const bool adaptive = mPrimAdaptivity > 1e-7 || mSecAdaptivity > 1e-7;
4886 const bool invertSurfaceOrientation = (!volume_to_mesh_internal::isBoolValue<InputValueType>()
4891 const InputTreeType& inputTree = inputGrid.tree();
4893 BoolTreeType intersectionTree(
false), adaptivityMask(
false);
4895 if (mAdaptivityMaskTree && mAdaptivityMaskTree->type() == BoolTreeType::treeType()) {
4896 const BoolTreeType *refAdaptivityMask=
4897 static_cast<const BoolTreeType*
>(mAdaptivityMaskTree.get());
4898 adaptivityMask.topologyUnion(*refAdaptivityMask);
4901 Int16TreeType signFlagsTree(0);
4907 volume_to_mesh_internal::identifySurfaceIntersectingVoxels(
4908 intersectionTree, inputTree, isovalue);
4910 volume_to_mesh_internal::applySurfaceMask(intersectionTree, adaptivityMask,
4911 inputGrid, mSurfaceMaskGrid, mInvertSurfaceMask, isovalue);
4913 if (intersectionTree.empty())
return;
4915 volume_to_mesh_internal::computeAuxiliaryData(
4916 signFlagsTree, pointIndexTree, intersectionTree, inputTree, isovalue);
4918 intersectionTree.clear();
4920 std::vector<Index32LeafNodeType*> pointIndexLeafNodes;
4921 pointIndexTree.getNodes(pointIndexLeafNodes);
4923 std::vector<Int16LeafNodeType*> signFlagsLeafNodes;
4924 signFlagsTree.getNodes(signFlagsLeafNodes);
4926 const tbb::blocked_range<size_t> auxiliaryLeafNodeRange(0, signFlagsLeafNodes.size());
4931 Int16TreeType* refSignFlagsTree =
nullptr;
4932 Index32TreeType* refPointIndexTree =
nullptr;
4933 InputTreeType
const* refInputTree =
nullptr;
4935 if (mRefGrid && mRefGrid->type() == InputGridType::gridType()) {
4937 const InputGridType* refGrid =
static_cast<const InputGridType*
>(mRefGrid.get());
4938 refInputTree = &refGrid->tree();
4940 if (!mRefSignTree && !mRefIdxTree) {
4944 typename Int16TreeType::Ptr refSignFlagsTreePt(
new Int16TreeType(0));
4945 typename Index32TreeType::Ptr refPointIndexTreePt(
4948 BoolTreeType refIntersectionTree(
false);
4950 volume_to_mesh_internal::identifySurfaceIntersectingVoxels(
4951 refIntersectionTree, *refInputTree, isovalue);
4953 volume_to_mesh_internal::computeAuxiliaryData(*refSignFlagsTreePt,
4954 *refPointIndexTreePt, refIntersectionTree, *refInputTree, isovalue);
4956 mRefSignTree = refSignFlagsTreePt;
4957 mRefIdxTree = refPointIndexTreePt;
4960 if (mRefSignTree && mRefIdxTree) {
4964 refSignFlagsTree =
static_cast<Int16TreeType*
>(mRefSignTree.get());
4965 refPointIndexTree =
static_cast<Index32TreeType*
>(mRefIdxTree.get());
4969 if (refSignFlagsTree && refPointIndexTree) {
4973 volume_to_mesh_internal::markSeamLineData(signFlagsTree, *refSignFlagsTree);
4975 if (mSeamPointListSize == 0) {
4979 std::vector<Int16LeafNodeType*> refSignFlagsLeafNodes;
4980 refSignFlagsTree->getNodes(refSignFlagsLeafNodes);
4982 std::unique_ptr<Index32[]> leafNodeOffsets(
4983 new Index32[refSignFlagsLeafNodes.size()]);
4986 volume_to_mesh_internal::LeafNodePointCount<Int16LeafNodeType::LOG2DIM>(
4987 refSignFlagsLeafNodes, leafNodeOffsets));
4991 for (
size_t n = 0, N = refSignFlagsLeafNodes.size(); n <
N; ++
n) {
4992 const Index32 tmp = leafNodeOffsets[
n];
4993 leafNodeOffsets[
n] =
count;
4996 mSeamPointListSize = size_t(count);
4999 if (mSeamPointListSize != 0) {
5001 mQuantizedSeamPoints.reset(
new uint32_t[mSeamPointListSize]);
5003 memset(mQuantizedSeamPoints.get(), 0,
sizeof(uint32_t) * mSeamPointListSize);
5005 std::vector<Index32LeafNodeType*> refPointIndexLeafNodes;
5006 refPointIndexTree->getNodes(refPointIndexLeafNodes);
5009 volume_to_mesh_internal::MapPoints<Index32LeafNodeType>(
5010 refPointIndexLeafNodes, refSignFlagsLeafNodes, leafNodeOffsets));
5014 if (mSeamPointListSize != 0) {
5017 volume_to_mesh_internal::SeamLineWeights<InputTreeType>(
5018 signFlagsLeafNodes, inputTree, *refPointIndexTree, *refSignFlagsTree,
5019 mQuantizedSeamPoints.get(), isovalue));
5024 const bool referenceMeshing = refSignFlagsTree && refPointIndexTree && refInputTree;
5029 std::unique_ptr<Index32[]> leafNodeOffsets(
new Index32[signFlagsLeafNodes.size()]);
5032 volume_to_mesh_internal::MergeVoxelRegions<InputGridType> mergeOp(
5033 inputGrid, pointIndexTree, pointIndexLeafNodes, signFlagsLeafNodes,
5034 isovalue, adaptivityThreshold, invertSurfaceOrientation);
5036 if (mAdaptivityGrid && mAdaptivityGrid->type() == FloatGridType::gridType()) {
5037 const FloatGridType* adaptivityGrid =
5038 static_cast<const FloatGridType*
>(mAdaptivityGrid.get());
5039 mergeOp.setSpatialAdaptivity(*adaptivityGrid);
5042 if (!adaptivityMask.empty()) {
5043 mergeOp.setAdaptivityMask(adaptivityMask);
5046 if (referenceMeshing) {
5047 mergeOp.setRefSignFlagsData(*refSignFlagsTree,
float(mSecAdaptivity));
5052 volume_to_mesh_internal::AdaptiveLeafNodePointCount<Index32LeafNodeType>
5053 op(pointIndexLeafNodes, signFlagsLeafNodes, leafNodeOffsets);
5059 volume_to_mesh_internal::LeafNodePointCount<Int16LeafNodeType::LOG2DIM>
5060 op(signFlagsLeafNodes, leafNodeOffsets);
5068 for (
size_t n = 0, N = signFlagsLeafNodes.size(); n <
N; ++
n) {
5069 const Index32 tmp = leafNodeOffsets[
n];
5074 mPointListSize = size_t(pointCount);
5076 mPointFlags.clear();
5083 volume_to_mesh_internal::ComputePoints<InputTreeType>
5084 op(mPoints.get(), inputTree, pointIndexLeafNodes,
5085 signFlagsLeafNodes, leafNodeOffsets,
transform, mIsovalue);
5087 if (referenceMeshing) {
5088 mPointFlags.resize(mPointListSize);
5089 op.setRefData(*refInputTree, *refPointIndexTree, *refSignFlagsTree,
5090 mQuantizedSeamPoints.get(), mPointFlags.data());
5099 mPolygonPoolListSize = signFlagsLeafNodes.size();
5100 mPolygons.reset(
new PolygonPool[mPolygonPoolListSize]);
5104 using PrimBuilder = volume_to_mesh_internal::AdaptivePrimBuilder;
5106 volume_to_mesh_internal::ComputePolygons<Int16TreeType, PrimBuilder>
5107 op(signFlagsLeafNodes, signFlagsTree, pointIndexTree,
5108 mPolygons, invertSurfaceOrientation);
5110 if (referenceMeshing) {
5111 op.setRefSignTree(refSignFlagsTree);
5118 using PrimBuilder = volume_to_mesh_internal::UniformPrimBuilder;
5120 volume_to_mesh_internal::ComputePolygons<Int16TreeType, PrimBuilder>
5121 op(signFlagsLeafNodes, signFlagsTree, pointIndexTree,
5122 mPolygons, invertSurfaceOrientation);
5124 if (referenceMeshing) {
5125 op.setRefSignTree(refSignFlagsTree);
5132 signFlagsTree.clear();
5133 pointIndexTree.clear();
5136 if (adaptive && mRelaxDisorientedTriangles) {
5137 volume_to_mesh_internal::relaxDisorientedTriangles(invertSurfaceOrientation,
5138 inputTree, transform, mPolygons, mPolygonPoolListSize, mPoints, mPointListSize);
5142 if (referenceMeshing) {
5143 volume_to_mesh_internal::subdivideNonplanarSeamLineQuads(
5144 mPolygons, mPolygonPoolListSize, mPoints, mPointListSize, mPointFlags);
5146 volume_to_mesh_internal::reviseSeamLineFlags(mPolygons, mPolygonPoolListSize, mPointFlags);
5159 template<
typename Gr
idType>
5162 const GridType& grid,
5163 std::vector<Vec3s>& points,
5164 std::vector<Vec3I>& triangles,
5165 std::vector<Vec4I>&
quads,
5168 bool relaxDisorientedTriangles)
5170 VolumeToMesh mesher(isovalue, adaptivity, relaxDisorientedTriangles);
5175 points.resize(mesher.pointListSize());
5178 volume_to_mesh_internal::PointListCopy ptnCpy(mesher.pointList(),
points);
5180 mesher.pointList().reset(
nullptr);
5186 size_t numQuads = 0, numTriangles = 0;
5187 for (
size_t n = 0, N = mesher.polygonPoolListSize(); n <
N; ++
n) {
5188 openvdb::tools::PolygonPool& polygons = polygonPoolList[
n];
5189 numTriangles += polygons.numTriangles();
5190 numQuads += polygons.numQuads();
5194 triangles.resize(numTriangles);
5196 quads.resize(numQuads);
5200 size_t qIdx = 0, tIdx = 0;
5201 for (
size_t n = 0, N = mesher.polygonPoolListSize(); n <
N; ++
n) {
5202 openvdb::tools::PolygonPool& polygons = polygonPoolList[
n];
5204 for (
size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
5205 quads[qIdx++] = polygons.quad(i);
5208 for (
size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
5209 triangles[tIdx++] = polygons.triangle(i);
5215 template<
typename Gr
idType>
5219 std::vector<Vec3s>&,
5220 std::vector<Vec3I>&,
5221 std::vector<Vec4I>&,
5226 OPENVDB_THROW(TypeError,
"volume to mesh conversion is supported only for scalar grids");
5233 template<
typename Gr
idType>
5236 const GridType& grid,
5237 std::vector<Vec3s>& points,
5238 std::vector<Vec3I>& triangles,
5239 std::vector<Vec4I>& quads,
5242 bool relaxDisorientedTriangles)
5244 doVolumeToMesh(grid, points, triangles, quads, isovalue, adaptivity, relaxDisorientedTriangles);
5248 template<
typename Gr
idType>
5251 const GridType& grid,
5252 std::vector<Vec3s>& points,
5253 std::vector<Vec4I>& quads,
5256 std::vector<Vec3I> triangles;
5257 doVolumeToMesh(grid, points, triangles, quads, isovalue, 0.0,
true);
5264 #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)
SharedPtr< TreeBase > Ptr
GLuint GLsizei const GLuint const GLintptr * offsets
**And then you can **find out if it s done
void reverse(I begin, I end)
#define OPENVDB_USE_VERSION_NAMESPACE
T dot(const Vec3< T > &v) const
Dot product.
GLfloat GLfloat GLfloat GLfloat GLfloat maxY
GLuint GLuint GLfloat weight
GLfloat GLfloat GLfloat GLfloat maxX
GLenum GLsizei GLsizei GLint * values
SIM_DerVector3 normalize() const
GLubyte GLubyte GLubyte GLubyte w
GLuint GLenum GLenum transform
Abstract base class for typed grids.
Mat3 transpose() const
returns transpose of this
SYS_FORCE_INLINE const_iterator end() const
std::vector< Index > IndexArray
void OIIO_API split(string_view str, std::vector< string_view > &result, string_view sep=string_view(), int maxsplit=-1)
GLint GLint GLsizei GLint GLenum GLenum type
ImageBuf OIIO_API min(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
IMATH_INTERNAL_NAMESPACE_HEADER_ENTER T abs(T a)
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) ...
GLfloat GLfloat GLfloat v2
ImageBuf OIIO_API max(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
GLdouble GLdouble GLdouble z
typedef int(WINAPI *PFNWGLRELEASEPBUFFERDCARBPROC)(HPBUFFERARB hPbuffer
OPENVDB_API const Index32 INVALID_IDX
Index64 pointCount(const PointDataTreeT &tree, const FilterT &filter=NullFilter(), const bool inCoreOnly=false, const bool threaded=true)
Count the total number of points in a PointDataTree.
Container class that associates a tree with a transform and metadata.
GLuint GLdouble GLdouble GLint GLint const GLdouble * points
Base class for typed trees.
GLuint GLsizei GLsizei * length
GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat maxZ
SharedPtr< const GridBase > ConstPtr
math::Vec3< Index32 > Vec3I
GA_API const UT_StringHolder N
SharedPtr< const TreeBase > ConstPtr
GLsizei const GLfloat * value
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
SIM_DerVector3 cross(const SIM_DerVector3 &lhs, const SIM_DerVector3 &rhs)
GLenum GLuint GLint GLenum face
#define OPENVDB_THROW(exception, message)