31 #ifndef OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
32 #define OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
34 #include <openvdb/Platform.h>
35 #include <openvdb/tree/ValueAccessor.h>
36 #include <openvdb/util/Util.h>
37 #include <openvdb/math/Operators.h>
38 #include <openvdb/tools/Morphology.h>
39 #include <openvdb/tree/LeafManager.h>
42 #include <boost/scoped_array.hpp>
43 #include <boost/scoped_ptr.hpp>
44 #include <boost/type_traits/is_scalar.hpp>
45 #include <boost/utility/enable_if.hpp>
46 #include <tbb/blocked_range.h>
47 #include <tbb/parallel_for.h>
48 #include <tbb/parallel_reduce.h>
77 template<
typename Gr
idType>
81 std::vector<Vec3s>& points,
82 std::vector<Vec4I>& quads,
83 double isovalue = 0.0);
96 template<
typename Gr
idType>
100 std::vector<Vec3s>& points,
101 std::vector<Vec3I>& triangles,
102 std::vector<Vec4I>& quads,
103 double isovalue = 0.0,
104 double adaptivity = 0.0);
120 inline PolygonPool(
const size_t numQuads,
const size_t numTriangles);
125 inline void resetQuads(
size_t size);
126 inline void clearQuads();
128 inline void resetTriangles(
size_t size);
129 inline void clearTriangles();
134 const size_t&
numQuads()
const {
return mNumQuads; }
149 const char&
quadFlags(
size_t n)
const {
return mQuadFlags[n]; }
158 inline bool trimQuads(
const size_t n,
bool reallocate =
false);
159 inline bool trimTrinagles(
const size_t n,
bool reallocate =
false);
165 size_t mNumQuads, mNumTriangles;
166 boost::scoped_array<openvdb::Vec4I> mQuads;
167 boost::scoped_array<openvdb::Vec3I> mTriangles;
168 boost::scoped_array<char> mQuadFlags, mTriangleFlags;
189 VolumeToMesh(
double isovalue = 0,
double adaptivity = 0);
196 const size_t& pointListSize()
const;
197 PointList& pointList();
199 const size_t& polygonPoolListSize()
const;
200 PolygonPoolList& polygonPoolList();
201 const PolygonPoolList& polygonPoolList()
const;
203 std::vector<unsigned char>& pointFlags();
204 const std::vector<unsigned char>& pointFlags()
const;
212 template<
typename Gr
idT>
213 void operator()(
const GridT&);
261 void partition(
unsigned partitions = 1,
unsigned activePart = 0);
266 PolygonPoolList mPolygons;
268 size_t mPointListSize, mSeamPointListSize, mPolygonPoolListSize;
269 double mIsovalue, mPrimAdaptivity, mSecAdaptivity;
276 bool mInvertSurfaceMask;
277 unsigned mPartitions, mActivePart;
279 boost::scoped_array<uint32_t> mQuantizedSeamPoints;
281 std::vector<unsigned char> mPointFlags;
295 const std::vector<Vec3d>& points,
296 const std::vector<Vec3d>& normals)
302 if (points.empty())
return avgPos;
304 for (
size_t n = 0, N = points.size(); n < N; ++n) {
308 avgPos /= double(points.size());
312 double m00=0,m01=0,m02=0,
319 for (
size_t n = 0, N = points.size(); n < N; ++n) {
321 const Vec3d& n_ref = normals[n];
324 m00 += n_ref[0] * n_ref[0];
325 m11 += n_ref[1] * n_ref[1];
326 m22 += n_ref[2] * n_ref[2];
328 m01 += n_ref[0] * n_ref[1];
329 m02 += n_ref[0] * n_ref[2];
330 m12 += n_ref[1] * n_ref[2];
333 rhs += n_ref * n_ref.
dot(points[n] - avgPos);
358 Mat3d D = Mat3d::identity();
361 double tolerance =
std::max(std::abs(eigenValues[0]), std::abs(eigenValues[1]));
362 tolerance =
std::max(tolerance, std::abs(eigenValues[2]));
366 for (
int i = 0; i < 3; ++i ) {
367 if (std::abs(eigenValues[i]) < tolerance) {
371 D[i][i] = 1.0 / eigenValues[i];
377 Mat3d pseudoInv = eigenVectors * D * eigenVectors.
transpose();
378 return avgPos + pseudoInv * rhs;
394 #ifdef OPENVDB_HAS_CXX11
395 typedef std::unique_ptr<T>
type;
409 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,
410 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,
411 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,
412 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,
413 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,
414 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,
415 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,
416 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};
421 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,
422 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,
423 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,
424 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,
425 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,
426 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,
427 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,
428 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};
435 {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},
436 {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},
437 {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},
438 {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},
439 {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},
440 {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},
441 {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},
442 {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},
443 {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},
444 {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},
445 {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},
446 {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},
447 {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},
448 {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},
449 {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},
450 {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},
451 {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},
452 {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},
453 {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},
454 {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},
455 {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},
456 {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},
457 {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},
458 {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},
459 {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},
460 {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},
461 {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},
462 {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},
463 {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},
464 {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},
465 {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},
466 {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},
467 {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},
468 {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},
469 {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},
470 {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},
471 {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},
472 {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},
473 {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},
474 {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},
475 {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},
476 {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},
477 {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},
478 {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},
479 {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},
480 {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},
481 {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},
482 {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},
483 {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},
484 {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},
485 {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},
486 {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},
487 {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},
488 {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},
489 {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},
490 {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},
491 {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},
492 {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},
493 {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},
494 {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},
495 {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},
496 {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},
497 {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},
498 {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},
499 {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},
500 {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},
501 {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},
502 {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},
503 {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},
504 {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},
505 {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},
506 {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},
507 {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},
508 {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},
509 {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},
510 {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},
511 {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},
512 {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},
513 {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},
514 {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},
515 {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},
516 {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},
517 {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},
518 {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},
519 {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},
520 {0,0,0,0,0,0,0,0,0,0,0,0,0}};
529 double epsilon = 0.001)
532 Vec3d normal = (p2-p0).cross(p1-p3);
534 const Vec3d centroid = (p0 + p1 + p2 + p3);
535 const double d = centroid.
dot(normal) * 0.25;
539 double absDist = std::abs(p0.dot(normal) - d);
540 if (absDist > epsilon)
return false;
542 absDist = std::abs(p1.dot(normal) - d);
543 if (absDist > epsilon)
return false;
545 absDist = std::abs(p2.dot(normal) - d);
546 if (absDist > epsilon)
return false;
548 absDist = std::abs(p3.dot(normal) - d);
549 if (absDist > epsilon)
return false;
573 assert(!(v.x() > 1.0) && !(v.y() > 1.0) && !(v.z() > 1.0));
574 assert(!(v.x() < 0.0) && !(v.y() < 0.0) && !(v.z() < 0.0));
589 v.y() = double(data & MASK_FIRST_10_BITS) * 0.0009775171;
591 v.x() = double(data & MASK_FIRST_10_BITS) * 0.0009775171;
603 template<
typename AccessorT>
605 evalCellSigns(
const AccessorT& accessor,
const Coord& ijk,
typename AccessorT::ValueType iso)
609 if (accessor.getValue(coord) < iso) signs |= 1u;
611 if (accessor.getValue(coord) < iso) signs |= 2u;
613 if (accessor.getValue(coord) < iso) signs |= 4u;
615 if (accessor.getValue(coord) < iso) signs |= 8u;
616 coord[1] += 1; coord[2] = ijk[2];
617 if (accessor.getValue(coord) < iso) signs |= 16u;
619 if (accessor.getValue(coord) < iso) signs |= 32u;
621 if (accessor.getValue(coord) < iso) signs |= 64u;
623 if (accessor.getValue(coord) < iso) signs |= 128u;
624 return uint8_t(signs);
630 template<
typename LeafT>
634 unsigned char signs = 0;
637 if (leaf.getValue(offset) < iso) signs |= 1u;
640 if (leaf.getValue(offset + 1) < iso) signs |= 8u;
643 if (leaf.getValue(offset + LeafT::DIM) < iso) signs |= 16u;
646 if (leaf.getValue(offset + LeafT::DIM + 1) < iso) signs |= 128u;
649 if (leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) ) < iso) signs |= 2u;
652 if (leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + 1) < iso) signs |= 4u;
655 if (leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM) < iso) signs |= 32u;
658 if (leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM + 1) < iso) signs |= 64u;
666 template<
class AccessorT>
669 const AccessorT& acc, Coord ijk,
typename AccessorT::ValueType iso)
673 if (sAmbiguousFace[
evalCellSigns(acc, ijk, iso)] == 3) signs = uint8_t(~signs);
674 }
else if (face == 3) {
676 if (sAmbiguousFace[
evalCellSigns(acc, ijk, iso)] == 1) signs = uint8_t(~signs);
677 }
else if (face == 2) {
679 if (sAmbiguousFace[
evalCellSigns(acc, ijk, iso)] == 4) signs = uint8_t(~signs);
680 }
else if (face == 4) {
682 if (sAmbiguousFace[
evalCellSigns(acc, ijk, iso)] == 2) signs = uint8_t(~signs);
683 }
else if (face == 5) {
685 if (sAmbiguousFace[
evalCellSigns(acc, ijk, iso)] == 6) signs = uint8_t(~signs);
686 }
else if (face == 6) {
688 if (sAmbiguousFace[
evalCellSigns(acc, ijk, iso)] == 5) signs = uint8_t(~signs);
693 template<
class AccessorT>
696 typename AccessorT::ValueType isovalue,
const int dim)
702 p[0] = accessor.getValue(coord) < isovalue;
704 p[1] = accessor.getValue(coord) < isovalue;
706 p[2] = accessor.getValue(coord) < isovalue;
708 p[3] = accessor.getValue(coord) < isovalue;
709 coord[1] += dim; coord[2] = ijk[2];
710 p[4] = accessor.getValue(coord) < isovalue;
712 p[5] = accessor.getValue(coord) < isovalue;
714 p[6] = accessor.getValue(coord) < isovalue;
716 p[7] = accessor.getValue(coord) < isovalue;
720 if (p[0]) signs |= 1u;
721 if (p[1]) signs |= 2u;
722 if (p[2]) signs |= 4u;
723 if (p[3]) signs |= 8u;
724 if (p[4]) signs |= 16u;
725 if (p[5]) signs |= 32u;
726 if (p[6]) signs |= 64u;
727 if (p[7]) signs |= 128u;
728 if (!sAdaptable[signs])
return true;
733 int i = ijk[0], ip = ijk[0] + hDim, ipp = ijk[0] + dim;
734 int j = ijk[1], jp = ijk[1] + hDim, jpp = ijk[1] + dim;
735 int k = ijk[2], kp = ijk[2] + hDim, kpp = ijk[2] + dim;
738 coord.reset(ip, j, k);
739 m = accessor.getValue(coord) < isovalue;
740 if (p[0] != m && p[1] != m)
return true;
743 coord.reset(ipp, j, kp);
744 m = accessor.getValue(coord) < isovalue;
745 if (p[1] != m && p[2] != m)
return true;
748 coord.reset(ip, j, kpp);
749 m = accessor.getValue(coord) < isovalue;
750 if (p[2] != m && p[3] != m)
return true;
753 coord.reset(i, j, kp);
754 m = accessor.getValue(coord) < isovalue;
755 if (p[0] != m && p[3] != m)
return true;
758 coord.reset(ip, jpp, k);
759 m = accessor.getValue(coord) < isovalue;
760 if (p[4] != m && p[5] != m)
return true;
763 coord.reset(ipp, jpp, kp);
764 m = accessor.getValue(coord) < isovalue;
765 if (p[5] != m && p[6] != m)
return true;
768 coord.reset(ip, jpp, kpp);
769 m = accessor.getValue(coord) < isovalue;
770 if (p[6] != m && p[7] != m)
return true;
773 coord.reset(i, jpp, kp);
774 m = accessor.getValue(coord) < isovalue;
775 if (p[7] != m && p[4] != m)
return true;
778 coord.reset(i, jp, k);
779 m = accessor.getValue(coord) < isovalue;
780 if (p[0] != m && p[4] != m)
return true;
783 coord.reset(ipp, jp, k);
784 m = accessor.getValue(coord) < isovalue;
785 if (p[1] != m && p[5] != m)
return true;
788 coord.reset(ipp, jp, kpp);
789 m = accessor.getValue(coord) < isovalue;
790 if (p[2] != m && p[6] != m)
return true;
794 coord.reset(i, jp, kpp);
795 m = accessor.getValue(coord) < isovalue;
796 if (p[3] != m && p[7] != m)
return true;
802 coord.reset(ip, jp, k);
803 m = accessor.getValue(coord) < isovalue;
804 if (p[0] != m && p[1] != m && p[4] != m && p[5] != m)
return true;
807 coord.reset(ipp, jp, kp);
808 m = accessor.getValue(coord) < isovalue;
809 if (p[1] != m && p[2] != m && p[5] != m && p[6] != m)
return true;
812 coord.reset(ip, jp, kpp);
813 m = accessor.getValue(coord) < isovalue;
814 if (p[2] != m && p[3] != m && p[6] != m && p[7] != m)
return true;
817 coord.reset(i, jp, kp);
818 m = accessor.getValue(coord) < isovalue;
819 if (p[0] != m && p[3] != m && p[4] != m && p[7] != m)
return true;
822 coord.reset(ip, j, kp);
823 m = accessor.getValue(coord) < isovalue;
824 if (p[0] != m && p[1] != m && p[2] != m && p[3] != m)
return true;
827 coord.reset(ip, jpp, kp);
828 m = accessor.getValue(coord) < isovalue;
829 if (p[4] != m && p[5] != m && p[6] != m && p[7] != m)
return true;
832 coord.reset(ip, jp, kp);
833 m = accessor.getValue(coord) < isovalue;
834 if (p[0] != m && p[1] != m && p[2] != m && p[3] != m &&
835 p[4] != m && p[5] != m && p[6] != m && p[7] != m)
return true;
844 template <
class LeafType>
846 mergeVoxels(LeafType& leaf,
const Coord& start,
int dim,
int regionId)
848 Coord ijk, end = start;
853 for (ijk[0] = start[0]; ijk[0] < end[0]; ++ijk[0]) {
854 for (ijk[1] = start[1]; ijk[1] < end[1]; ++ijk[1]) {
855 for (ijk[2] = start[2]; ijk[2] < end[2]; ++ijk[2]) {
856 leaf.setValueOnly(ijk, regionId);
865 template <
class LeafType>
868 typename LeafType::ValueType::value_type adaptivity)
870 if (adaptivity < 1e-6)
return false;
872 typedef typename LeafType::ValueType VecT;
873 Coord ijk, end = start;
878 std::vector<VecT> norms;
879 for (ijk[0] = start[0]; ijk[0] < end[0]; ++ijk[0]) {
880 for (ijk[1] = start[1]; ijk[1] < end[1]; ++ijk[1]) {
881 for (ijk[2] = start[2]; ijk[2] < end[2]; ++ijk[2]) {
883 if(!leaf.isValueOn(ijk))
continue;
884 norms.push_back(leaf.getValue(ijk));
889 size_t N = norms.size();
890 for (
size_t ni = 0; ni < N; ++ni) {
891 VecT n_i = norms[ni];
892 for (
size_t nj = 0; nj < N; ++nj) {
893 VecT n_j = norms[nj];
894 if ((1.0 - n_i.dot(n_j)) > adaptivity)
return false;
904 template<
typename TreeT,
typename LeafManagerT>
908 typedef typename TreeT::ValueType
ValueT;
911 typedef typename TreeT::template ValueConverter<int>::Type
IntTreeT;
914 typedef typename TreeT::template ValueConverter<Int16>::Type
Int16TreeT;
920 SignData(
const TreeT& distTree,
const LeafManagerT& leafs, ValueT iso);
922 void run(
bool threaded =
true);
924 typename Int16TreeT::Ptr
signTree()
const {
return mSignTree; }
925 typename IntTreeT::Ptr
idxTree()
const {
return mIdxTree; }
930 void operator()(
const tbb::blocked_range<size_t>&);
933 mSignTree->merge(*rhs.mSignTree);
934 mIdxTree->merge(*rhs.mIdxTree);
939 const TreeT& mDistTree;
942 const LeafManagerT& mLeafs;
945 typename Int16TreeT::Ptr mSignTree;
946 Int16AccessorT mSignAcc;
948 typename IntTreeT::Ptr mIdxTree;
949 IntAccessorT mIdxAcc;
954 template<
typename TreeT,
typename LeafManagerT>
956 const LeafManagerT& leafs,
ValueT iso)
957 : mDistTree(distTree)
958 , mDistAcc(mDistTree)
962 , mSignAcc(*mSignTree)
969 template<
typename TreeT,
typename LeafManagerT>
971 : mDistTree(rhs.mDistTree)
972 , mDistAcc(mDistTree)
974 , mIsovalue(rhs.mIsovalue)
976 , mSignAcc(*mSignTree)
983 template<
typename TreeT,
typename LeafManagerT>
987 if (threaded) tbb::parallel_reduce(mLeafs.getRange(), *
this);
988 else (*
this)(mLeafs.getRange());
991 template<
typename TreeT,
typename LeafManagerT>
995 typedef typename Int16TreeT::LeafNodeType Int16LeafT;
996 typedef typename IntTreeT::LeafNodeType IntLeafT;
997 typename LeafManagerT::TreeType::LeafNodeType::ValueOnCIter iter;
998 unsigned char signs, face;
1003 for (
size_t n = range.begin(); n != range.end(); ++n) {
1005 bool collectedData =
false;
1007 coord = mLeafs.leaf(n).origin();
1009 if (!signLeafPt.get()) signLeafPt.reset(
new Int16LeafT(coord, 0));
1010 else signLeafPt->setOrigin(coord);
1012 const typename TreeT::LeafNodeType *leafPt = mDistAcc.probeConstLeaf(coord);
1014 coord.offset(TreeT::LeafNodeType::DIM - 1);
1016 for (iter = mLeafs.leaf(n).cbeginValueOn(); iter; ++iter) {
1018 ijk = iter.getCoord();
1020 if (leafPt && ijk[0] < coord[0] && ijk[1] < coord[1] && ijk[2] < coord[2]) {
1026 if (signs != 0 && signs != 0xFF) {
1029 if (
bool(signs & 0x1) != bool(signs & 0x2)) flags |=
XEDGE;
1030 if (
bool(signs & 0x1) != bool(signs & 0x10)) flags |=
YEDGE;
1031 if (
bool(signs & 0x1) != bool(signs & 0x8)) flags |=
ZEDGE;
1038 signLeafPt->setValue(ijk, flags);
1039 collectedData =
true;
1043 if (collectedData) {
1045 IntLeafT* idxLeaf = mIdxAcc.touchLeaf(coord);
1046 idxLeaf->topologyUnion(*signLeafPt);
1047 typename IntLeafT::ValueOnIter it = idxLeaf->beginValueOn();
1052 mSignAcc.addLeaf(signLeafPt.release());
1065 CountPoints(std::vector<size_t>& pointList) : mPointList(pointList) {}
1067 template <
typename LeafNodeType>
1072 typename LeafNodeType::ValueOnCIter iter = leaf.cbeginValueOn();
1073 for (; iter; ++iter) {
1074 points += size_t(sEdgeGroupTable[(
SIGNS & iter.getValue())][0]);
1077 mPointList[leafIndex] = points;
1081 std::vector<size_t>& mPointList;
1086 template<
typename Int16TreeT>
1092 MapPoints(std::vector<size_t>& pointList,
const Int16TreeT& signTree)
1093 : mPointList(pointList)
1094 , mSignAcc(signTree)
1098 template <
typename LeafNodeType>
1101 size_t ptnIdx = mPointList[leafIndex];
1102 typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
1104 const typename Int16TreeT::LeafNodeType *signLeafPt =
1107 for (; iter; ++iter) {
1108 iter.setValue(static_cast<typename LeafNodeType::ValueType>(ptnIdx));
1109 unsigned signs =
SIGNS & signLeafPt->getValue(iter.pos());
1110 ptnIdx += size_t(sEdgeGroupTable[signs][0]);
1115 std::vector<size_t>& mPointList;
1116 Int16AccessorT mSignAcc;
1121 template<
typename IntTreeT>
1134 template <
typename LeafNodeType>
1142 typename IntLeafT::ValueOnIter iter = tmpLeaf.beginValueOn();
1143 for (; iter; ++iter) {
1144 if(iter.getValue() == 0) {
1146 regions += size_t(sEdgeGroupTable[(
SIGNS & leaf.getValue(iter.pos()))][0]);
1150 int onVoxelCount = int(tmpLeaf.onVoxelCount());
1151 while (onVoxelCount > 0) {
1153 iter = tmpLeaf.beginValueOn();
1154 int regionId = iter.getValue();
1155 for (; iter; ++iter) {
1156 if (iter.getValue() == regionId) {
1163 mRegions[leafIndex] = regions;
1167 IntAccessorT mIdxAcc;
1168 std::vector<size_t>& mRegions;
1176 inline double evalRoot(
double v0,
double v1,
double iso) {
return (iso - v0) / (v1 - v0); }
1180 template<
typename LeafT>
1184 values[0] = double(leaf.getValue(offset));
1185 values[3] = double(leaf.getValue(offset + 1));
1186 values[4] = double(leaf.getValue(offset + LeafT::DIM));
1187 values[7] = double(leaf.getValue(offset + LeafT::DIM + 1));
1188 values[1] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM)));
1189 values[2] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + 1));
1190 values[5] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM));
1191 values[6] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM + 1));
1196 template<
typename AccessorT>
1201 values[0] = double(acc.getValue(coord));
1204 values[1] = double(acc.getValue(coord));
1207 values[2] = double(acc.getValue(coord));
1210 values[3] = double(acc.getValue(coord));
1212 coord[1] += 1; coord[2] = ijk[2];
1213 values[4] = double(acc.getValue(coord));
1216 values[5] = double(acc.getValue(coord));
1219 values[6] = double(acc.getValue(coord));
1222 values[7] = double(acc.getValue(coord));
1229 unsigned char edgeGroup,
double iso)
1231 Vec3d avg(0.0, 0.0, 0.0);
1234 if (sEdgeGroupTable[signs][1] == edgeGroup) {
1235 avg[0] +=
evalRoot(values[0], values[1], iso);
1239 if (sEdgeGroupTable[signs][2] == edgeGroup) {
1241 avg[2] +=
evalRoot(values[1], values[2], iso);
1245 if (sEdgeGroupTable[signs][3] == edgeGroup) {
1246 avg[0] +=
evalRoot(values[3], values[2], iso);
1251 if (sEdgeGroupTable[signs][4] == edgeGroup) {
1252 avg[2] +=
evalRoot(values[0], values[3], iso);
1256 if (sEdgeGroupTable[signs][5] == edgeGroup) {
1257 avg[0] +=
evalRoot(values[4], values[5], iso);
1262 if (sEdgeGroupTable[signs][6] == edgeGroup) {
1265 avg[2] +=
evalRoot(values[5], values[6], iso);
1269 if (sEdgeGroupTable[signs][7] == edgeGroup) {
1270 avg[0] +=
evalRoot(values[7], values[6], iso);
1276 if (sEdgeGroupTable[signs][8] == edgeGroup) {
1278 avg[2] +=
evalRoot(values[4], values[7], iso);
1282 if (sEdgeGroupTable[signs][9] == edgeGroup) {
1283 avg[1] +=
evalRoot(values[0], values[4], iso);
1287 if (sEdgeGroupTable[signs][10] == edgeGroup) {
1289 avg[1] +=
evalRoot(values[1], values[5], iso);
1293 if (sEdgeGroupTable[signs][11] == edgeGroup) {
1295 avg[1] +=
evalRoot(values[2], values[6], iso);
1300 if (sEdgeGroupTable[signs][12] == edgeGroup) {
1301 avg[1] +=
evalRoot(values[3], values[7], iso);
1307 double w = 1.0 / double(samples);
1321 unsigned char signsMask,
unsigned char edgeGroup,
double iso)
1323 avg =
Vec3d(0.0, 0.0, 0.0);
1326 if (sEdgeGroupTable[signs][1] == edgeGroup
1327 && sEdgeGroupTable[signsMask][1] == 0) {
1328 avg[0] +=
evalRoot(values[0], values[1], iso);
1332 if (sEdgeGroupTable[signs][2] == edgeGroup
1333 && sEdgeGroupTable[signsMask][2] == 0) {
1335 avg[2] +=
evalRoot(values[1], values[2], iso);
1339 if (sEdgeGroupTable[signs][3] == edgeGroup
1340 && sEdgeGroupTable[signsMask][3] == 0) {
1341 avg[0] +=
evalRoot(values[3], values[2], iso);
1346 if (sEdgeGroupTable[signs][4] == edgeGroup
1347 && sEdgeGroupTable[signsMask][4] == 0) {
1348 avg[2] +=
evalRoot(values[0], values[3], iso);
1352 if (sEdgeGroupTable[signs][5] == edgeGroup
1353 && sEdgeGroupTable[signsMask][5] == 0) {
1354 avg[0] +=
evalRoot(values[4], values[5], iso);
1359 if (sEdgeGroupTable[signs][6] == edgeGroup
1360 && sEdgeGroupTable[signsMask][6] == 0) {
1363 avg[2] +=
evalRoot(values[5], values[6], iso);
1367 if (sEdgeGroupTable[signs][7] == edgeGroup
1368 && sEdgeGroupTable[signsMask][7] == 0) {
1369 avg[0] +=
evalRoot(values[7], values[6], iso);
1375 if (sEdgeGroupTable[signs][8] == edgeGroup
1376 && sEdgeGroupTable[signsMask][8] == 0) {
1378 avg[2] +=
evalRoot(values[4], values[7], iso);
1382 if (sEdgeGroupTable[signs][9] == edgeGroup
1383 && sEdgeGroupTable[signsMask][9] == 0) {
1384 avg[1] +=
evalRoot(values[0], values[4], iso);
1388 if (sEdgeGroupTable[signs][10] == edgeGroup
1389 && sEdgeGroupTable[signsMask][10] == 0) {
1391 avg[1] +=
evalRoot(values[1], values[5], iso);
1395 if (sEdgeGroupTable[signs][11] == edgeGroup
1396 && sEdgeGroupTable[signsMask][11] == 0) {
1398 avg[1] +=
evalRoot(values[2], values[6], iso);
1403 if (sEdgeGroupTable[signs][12] == edgeGroup
1404 && sEdgeGroupTable[signsMask][12] == 0) {
1405 avg[1] +=
evalRoot(values[3], values[7], iso);
1411 double w = 1.0 / double(samples);
1425 unsigned char signs,
unsigned char edgeGroup,
double iso)
1427 std::vector<Vec3d> samples;
1430 std::vector<double> weights;
1433 Vec3d avg(0.0, 0.0, 0.0);
1435 if (sEdgeGroupTable[signs][1] == edgeGroup) {
1436 avg[0] =
evalRoot(values[0], values[1], iso);
1440 samples.push_back(avg);
1441 weights.push_back((avg-p).lengthSqr());
1444 if (sEdgeGroupTable[signs][2] == edgeGroup) {
1447 avg[2] =
evalRoot(values[1], values[2], iso);
1449 samples.push_back(avg);
1450 weights.push_back((avg-p).lengthSqr());
1453 if (sEdgeGroupTable[signs][3] == edgeGroup) {
1454 avg[0] =
evalRoot(values[3], values[2], iso);
1458 samples.push_back(avg);
1459 weights.push_back((avg-p).lengthSqr());
1462 if (sEdgeGroupTable[signs][4] == edgeGroup) {
1465 avg[2] =
evalRoot(values[0], values[3], iso);
1467 samples.push_back(avg);
1468 weights.push_back((avg-p).lengthSqr());
1471 if (sEdgeGroupTable[signs][5] == edgeGroup) {
1472 avg[0] =
evalRoot(values[4], values[5], iso);
1476 samples.push_back(avg);
1477 weights.push_back((avg-p).lengthSqr());
1480 if (sEdgeGroupTable[signs][6] == edgeGroup) {
1483 avg[2] =
evalRoot(values[5], values[6], iso);
1485 samples.push_back(avg);
1486 weights.push_back((avg-p).lengthSqr());
1489 if (sEdgeGroupTable[signs][7] == edgeGroup) {
1490 avg[0] =
evalRoot(values[7], values[6], iso);
1494 samples.push_back(avg);
1495 weights.push_back((avg-p).lengthSqr());
1498 if (sEdgeGroupTable[signs][8] == edgeGroup) {
1501 avg[2] =
evalRoot(values[4], values[7], iso);
1503 samples.push_back(avg);
1504 weights.push_back((avg-p).lengthSqr());
1507 if (sEdgeGroupTable[signs][9] == edgeGroup) {
1509 avg[1] =
evalRoot(values[0], values[4], iso);
1512 samples.push_back(avg);
1513 weights.push_back((avg-p).lengthSqr());
1516 if (sEdgeGroupTable[signs][10] == edgeGroup) {
1518 avg[1] =
evalRoot(values[1], values[5], iso);
1521 samples.push_back(avg);
1522 weights.push_back((avg-p).lengthSqr());
1525 if (sEdgeGroupTable[signs][11] == edgeGroup) {
1527 avg[1] =
evalRoot(values[2], values[6], iso);
1530 samples.push_back(avg);
1531 weights.push_back((avg-p).lengthSqr());
1534 if (sEdgeGroupTable[signs][12] == edgeGroup) {
1536 avg[1] =
evalRoot(values[3], values[7], iso);
1539 samples.push_back(avg);
1540 weights.push_back((avg-p).lengthSqr());
1547 for (
size_t i = 0, I = weights.size(); i < I; ++i) {
1548 minWeight =
std::min(minWeight, weights[i]);
1549 maxWeight =
std::max(maxWeight, weights[i]);
1552 const double offset = maxWeight + minWeight * 0.1;
1553 for (
size_t i = 0, I = weights.size(); i < I; ++i) {
1554 weights[i] = offset - weights[i];
1558 double weightSum = 0.0;
1559 for (
size_t i = 0, I = weights.size(); i < I; ++i) {
1560 weightSum += weights[i];
1567 if (samples.size() > 1) {
1568 for (
size_t i = 0, I = samples.size(); i < I; ++i) {
1569 avg += samples[i] * (weights[i] / weightSum);
1572 avg = samples.front();
1583 const std::vector<double>& values,
unsigned char signs,
double iso)
1585 for (
size_t n = 1, N = sEdgeGroupTable[signs][0] + 1; n < N; ++n) {
1586 points.push_back(
computePoint(values, signs, uint8_t(n), iso));
1595 matchEdgeGroup(
unsigned char groupId,
unsigned char lhsSigns,
unsigned char rhsSigns)
1598 for (
size_t i = 1; i <= 12; ++i) {
1599 if (sEdgeGroupTable[lhsSigns][i] == groupId && sEdgeGroupTable[rhsSigns][i] != 0) {
1600 id = sEdgeGroupTable[rhsSigns][i];
1614 const std::vector<double>& lhsValues,
const std::vector<double>& rhsValues,
1615 unsigned char lhsSigns,
unsigned char rhsSigns,
1616 double iso,
size_t pointIdx,
const boost::scoped_array<uint32_t>& seamPoints)
1618 for (
size_t n = 1, N = sEdgeGroupTable[lhsSigns][0] + 1; n < N; ++n) {
1624 const unsigned char e = uint8_t(
id);
1625 uint32_t& quantizedPoint = seamPoints[pointIdx + (
id - 1)];
1630 weightedPointMask.push_back(
true);
1632 points.push_back(
computePoint(rhsValues, rhsSigns, e, iso));
1633 weightedPointMask.push_back(
false);
1637 points.push_back(
computePoint(lhsValues, lhsSigns, uint8_t(n), iso));
1638 weightedPointMask.push_back(
false);
1644 template <
typename TreeT,
typename LeafManagerT>
1650 typedef typename TreeT::template ValueConverter<int>::Type
IntTreeT;
1654 typedef typename TreeT::template ValueConverter<Int16>::Type
Int16TreeT;
1662 GenPoints(
const LeafManagerT& signLeafs,
const TreeT& distTree,
1663 IntTreeT& idxTree, PointList& points, std::vector<size_t>& indices,
1666 void run(
bool threaded =
true);
1668 void setRefData(
const Int16TreeT* refSignTree = NULL,
const TreeT* refDistTree = NULL,
1669 IntTreeT* refIdxTree = NULL,
const QuantizedPointList* seamPoints = NULL,
1670 std::vector<unsigned char>* mSeamPointMaskPt = NULL);
1675 void operator()(
const tbb::blocked_range<size_t>&)
const;
1678 const LeafManagerT& mSignLeafs;
1684 std::vector<size_t>& mIndices;
1686 const double mIsovalue;
1689 const Int16TreeT *mRefSignTreePt;
1690 const TreeT* mRefDistTreePt;
1691 const IntTreeT* mRefIdxTreePt;
1692 const QuantizedPointList* mSeamPointsPt;
1693 std::vector<unsigned char>* mSeamPointMaskPt;
1697 template <
typename TreeT,
typename LeafManagerT>
1699 const TreeT& distTree,
IntTreeT& idxTree, PointList& points,
1700 std::vector<size_t>& indices,
const math::Transform& xform,
double iso)
1701 : mSignLeafs(signLeafs)
1702 , mDistAcc(distTree)
1708 , mRefSignTreePt(NULL)
1709 , mRefDistTreePt(NULL)
1710 , mRefIdxTreePt(NULL)
1711 , mSeamPointsPt(NULL)
1712 , mSeamPointMaskPt(NULL)
1717 template <
typename TreeT,
typename LeafManagerT>
1721 if (threaded) tbb::parallel_for(mSignLeafs.getRange(), *
this);
1722 else (*
this)(mSignLeafs.getRange());
1726 template <
typename TreeT,
typename LeafManagerT>
1730 const TreeT *refDistTree,
1733 std::vector<unsigned char>* seamPointMask)
1735 mRefSignTreePt = refSignTree;
1736 mRefDistTreePt = refDistTree;
1737 mRefIdxTreePt = refIdxTree;
1738 mSeamPointsPt = seamPoints;
1739 mSeamPointMaskPt = seamPointMask;
1743 template <
typename TreeT,
typename LeafManagerT>
1747 typename IntTreeT::LeafNodeType::ValueOnIter iter;
1748 unsigned char signs, refSigns;
1751 std::vector<Vec3d> points(4);
1752 std::vector<bool> weightedPointMask(4);
1753 std::vector<double> values(8), refValues(8);
1759 boost::scoped_ptr<Int16CAccessorT> refSignAcc;
1760 if (mRefSignTreePt) refSignAcc.reset(
new Int16CAccessorT(*mRefSignTreePt));
1762 boost::scoped_ptr<IntCAccessorT> refIdxAcc;
1763 if (mRefIdxTreePt) refIdxAcc.reset(
new IntCAccessorT(*mRefIdxTreePt));
1765 boost::scoped_ptr<AccessorT> refDistAcc;
1766 if (mRefDistTreePt) refDistAcc.reset(
new AccessorT(*mRefDistTreePt));
1769 for (
size_t n = range.begin(); n != range.end(); ++n) {
1771 coord = mSignLeafs.leaf(n).origin();
1773 const typename TreeT::LeafNodeType *leafPt = mDistAcc.probeConstLeaf(coord);
1774 typename IntTreeT::LeafNodeType *idxLeafPt = idxAcc.
probeLeaf(coord);
1778 const typename Int16TreeT::LeafNodeType *refSignLeafPt = NULL;
1779 if (refSignAcc) refSignLeafPt = refSignAcc->probeConstLeaf(coord);
1781 const typename IntTreeT::LeafNodeType *refIdxLeafPt = NULL;
1782 if (refIdxAcc) refIdxLeafPt = refIdxAcc->probeConstLeaf(coord);
1784 const typename TreeT::LeafNodeType *refDistLeafPt = NULL;
1785 if (refDistAcc) refDistLeafPt = refDistAcc->probeConstLeaf(coord);
1789 size_t ptnIdx = mIndices[n];
1790 coord.offset(TreeT::LeafNodeType::DIM - 1);
1794 for (iter = idxLeafPt->beginValueOn(); iter; ++iter) {
1796 if(iter.getValue() != 0)
continue;
1798 iter.setValue(static_cast<typename IntTreeT::ValueType>(ptnIdx));
1800 offset = iter.pos();
1801 ijk = iter.getCoord();
1803 const bool inclusiveCell = ijk[0] < coord[0] && ijk[1] < coord[1] && ijk[2] < coord[2];
1805 const Int16& flags = mSignLeafs.leaf(n).getValue(offset);
1806 signs = uint8_t(
SIGNS & flags);
1809 if ((flags &
SEAM) && refSignLeafPt && refIdxLeafPt) {
1810 if (refSignLeafPt->isValueOn(offset)) {
1811 refSigns = uint8_t(
SIGNS & refSignLeafPt->getValue(offset));
1821 weightedPointMask.clear();
1823 if (refSigns == 0) {
1830 computeCellPoints(points, weightedPointMask, values, refValues, signs, refSigns,
1831 mIsovalue, refIdxLeafPt->getValue(offset), *mSeamPointsPt);
1835 for (
size_t i = 0, I = points.size(); i < I; ++i) {
1838 points[i][0] += double(ijk[0]);
1839 points[i][1] += double(ijk[1]);
1840 points[i][2] += double(ijk[2]);
1843 points[i] = mTransform.indexToWorld(points[i]);
1845 mPoints[ptnIdx][0] = float(points[i][0]);
1846 mPoints[ptnIdx][1] = float(points[i][1]);
1847 mPoints[ptnIdx][2] = float(points[i][2]);
1849 if (mSeamPointMaskPt && !weightedPointMask.empty() && weightedPointMask[i]) {
1850 (*mSeamPointMaskPt)[ptnIdx] = 1;
1858 int onVoxelCount = int(idxLeafPt->onVoxelCount());
1859 while (onVoxelCount > 0) {
1861 iter = idxLeafPt->beginValueOn();
1862 int regionId = iter.getValue(), count = 0;
1864 Vec3d avg(0.0), point;
1866 for (; iter; ++iter) {
1867 if (iter.getValue() != regionId)
continue;
1869 iter.setValue(static_cast<typename IntTreeT::ValueType>(ptnIdx));
1873 ijk = iter.getCoord();
1874 offset = iter.pos();
1876 signs = uint8_t(
SIGNS & mSignLeafs.leaf(n).getValue(offset));
1878 if (ijk[0] < coord[0] && ijk[1] < coord[1] && ijk[2] < coord[2]) {
1887 avg[0] += double(ijk[0]) + points[0][0];
1888 avg[1] += double(ijk[1]) + points[0][1];
1889 avg[2] += double(ijk[2]) + points[0][2];
1896 double w = 1.0 / double(count);
1902 avg = mTransform.indexToWorld(avg);
1904 mPoints[ptnIdx][0] = float(avg[0]);
1905 mPoints[ptnIdx][1] = float(avg[1]);
1906 mPoints[ptnIdx][2] = float(avg[2]);
1917 template<
typename TreeT>
1923 typedef typename TreeT::template ValueConverter<int>::Type
IntTreeT;
1926 typedef typename TreeT::template ValueConverter<Int16>::Type
Int16TreeT;
1933 SeamWeights(
const TreeT& distTree,
const Int16TreeT& refSignTree,
1934 IntTreeT& refIdxTree, QuantizedPointList& points,
double iso);
1936 template <
typename LeafNodeType>
1937 void operator()(LeafNodeType &signLeaf,
size_t leafIndex)
const;
1941 Int16AccessorT mRefSignAcc;
1942 IntAccessorT mRefIdxAcc;
1944 QuantizedPointList& mPoints;
1945 const double mIsovalue;
1949 template<
typename TreeT>
1952 : mDistAcc(distTree)
1953 , mRefSignAcc(refSignTree)
1954 , mRefIdxAcc(refIdxTree)
1961 template<
typename TreeT>
1962 template <
typename LeafNodeType>
1966 Coord coord = signLeaf.origin();
1967 const typename Int16TreeT::LeafNodeType *refSignLeafPt = mRefSignAcc.probeConstLeaf(coord);
1969 if (!refSignLeafPt)
return;
1971 const typename TreeT::LeafNodeType *distLeafPt = mDistAcc.probeConstLeaf(coord);
1972 const typename IntTreeT::LeafNodeType *refIdxLeafPt = mRefIdxAcc.probeConstLeaf(coord);
1974 std::vector<double> values(8);
1975 unsigned char lhsSigns, rhsSigns;
1980 coord.offset(TreeT::LeafNodeType::DIM - 1);
1982 typename LeafNodeType::ValueOnCIter iter = signLeaf.cbeginValueOn();
1983 for (; iter; ++iter) {
1985 offset = iter.pos();
1986 ijk = iter.getCoord();
1988 const bool inclusiveCell = ijk[0] < coord[0] && ijk[1] < coord[1] && ijk[2] < coord[2];
1990 if ((iter.getValue() &
SEAM) && refSignLeafPt->isValueOn(offset)) {
1992 lhsSigns = uint8_t(
SIGNS & iter.getValue());
1993 rhsSigns = uint8_t(
SIGNS & refSignLeafPt->getValue(offset));
1996 if (inclusiveCell) {
2003 for (
size_t n = 1, N = sEdgeGroupTable[lhsSigns][0] + 1; n < N; ++n) {
2009 uint32_t& data = mPoints[refIdxLeafPt->getValue(offset) + (
id - 1)];
2014 point, values, lhsSigns, rhsSigns, uint8_t(n), mIsovalue);
2016 if (smaples > 0) data =
packPoint(point);
2031 template <
typename TreeT,
typename LeafManagerT>
2038 typedef typename TreeT::template ValueConverter<int>::Type
IntTreeT;
2041 typedef typename TreeT::template ValueConverter<bool>::Type
BoolTreeT;
2043 typedef typename LeafManagerT::TreeType::template ValueConverter<Int16>::Type
Int16TreeT;
2046 typedef typename TreeT::template ValueConverter<float>::Type
FloatTreeT;
2053 const TreeT& distTree, IntTreeT& idxTree, ValueT iso, ValueT adaptivity);
2055 void run(
bool threaded =
true);
2058 const math::Transform& distGridXForm,
const FloatGridT& adaptivityField);
2062 void setRefData(
const Int16TreeT* signTree, ValueT adaptivity);
2067 void operator()(
const tbb::blocked_range<size_t>&)
const;
2071 const LeafManagerT& mSignLeafs;
2073 const Int16TreeT& mSignTree;
2074 Int16AccessorT mSignAcc;
2076 const TreeT& mDistTree;
2080 ValueT mIsovalue, mSurfaceAdaptivity, mInternalAdaptivity;
2083 const FloatGridT* mAdaptivityGrid;
2084 const BoolTreeT* mAdaptivityMask;
2086 const Int16TreeT* mRefSignTree;
2090 template <
typename TreeT,
typename LeafManagerT>
2092 const LeafManagerT& signLeafs,
const Int16TreeT& signTree,
2094 : mSignLeafs(signLeafs)
2095 , mSignTree(signTree)
2096 , mSignAcc(mSignTree)
2097 , mDistTree(distTree)
2098 , mDistAcc(mDistTree)
2101 , mSurfaceAdaptivity(adaptivity)
2102 , mInternalAdaptivity(adaptivity)
2104 , mAdaptivityGrid(NULL)
2105 , mAdaptivityMask(NULL)
2106 , mRefSignTree(NULL)
2111 template <
typename TreeT,
typename LeafManagerT>
2115 if (threaded) tbb::parallel_for(mSignLeafs.getRange(), *
this);
2116 else (*
this)(mSignLeafs.getRange());
2120 template <
typename TreeT,
typename LeafManagerT>
2125 mTransform = &distGridXForm;
2126 mAdaptivityGrid = &adaptivityField;
2130 template <
typename TreeT,
typename LeafManagerT>
2134 mAdaptivityMask = mask;
2137 template <
typename TreeT,
typename LeafManagerT>
2141 mRefSignTree = signTree;
2142 mInternalAdaptivity = adaptivity;
2146 template <
typename TreeT,
typename LeafManagerT>
2152 typedef typename TreeT::LeafNodeType LeafT;
2153 typedef typename IntTreeT::LeafNodeType IntLeafT;
2154 typedef typename BoolTreeT::LeafNodeType BoolLeafT;
2155 typedef typename LeafT::template ValueConverter<Vec3T>::Type Vec3LeafT;
2157 const int LeafDim = LeafT::DIM;
2161 typename LeafManagerT::TreeType::LeafNodeType::ValueOnCIter iter;
2164 boost::scoped_ptr<FloatTreeCAccessorT> adaptivityAcc;
2165 if (mAdaptivityGrid) {
2166 adaptivityAcc.reset(
new FloatTreeCAccessorT(mAdaptivityGrid->tree()));
2170 boost::scoped_ptr<Int16TreeCAccessorT> refAcc;
2172 refAcc.reset(
new Int16TreeCAccessorT(*mRefSignTree));
2176 boost::scoped_ptr<BoolTreeCAccessorT> maskAcc;
2177 if (mAdaptivityMask) {
2178 maskAcc.reset(
new BoolTreeCAccessorT(*mAdaptivityMask));
2183 Vec3LeafT gradients;
2186 for (
size_t n = range.begin(); n != range.end(); ++n) {
2188 mask.setValuesOff();
2190 const Coord& origin = mSignLeafs.leaf(n).origin();
2192 ValueT adaptivity = (refAcc && !refAcc->probeConstLeaf(origin)) ?
2193 mInternalAdaptivity : mSurfaceAdaptivity;
2195 IntLeafT& idxLeaf = *idxAcc.
probeLeaf(origin);
2197 end[0] = origin[0] + LeafDim;
2198 end[1] = origin[1] + LeafDim;
2199 end[2] = origin[2] + LeafDim;
2203 const BoolLeafT* maskLeaf = maskAcc->probeConstLeaf(origin);
2204 if (maskLeaf != NULL) {
2205 typename BoolLeafT::ValueOnCIter it;
2206 for (it = maskLeaf->cbeginValueOn(); it; ++it) {
2207 mask.setActiveState(it.getCoord() & ~1u,
true);
2213 LeafT adaptivityLeaf(origin, adaptivity);
2214 if (mAdaptivityGrid) {
2215 for (
Index offset = 0; offset < LeafT::NUM_VALUES; ++offset) {
2216 ijk = adaptivityLeaf.offsetToGlobalCoord(offset);
2217 Vec3d xyz = mAdaptivityGrid->transform().worldToIndex(
2218 mTransform->indexToWorld(ijk));
2220 adaptivityLeaf.setValueOnly(offset, tmpA * adaptivity);
2225 for (iter = mSignLeafs.leaf(n).cbeginValueOn(); iter; ++iter) {
2226 unsigned char signs =
static_cast<unsigned char>(
SIGNS & int(iter.getValue()));
2227 if (!sAdaptable[signs] || sEdgeGroupTable[signs][0] > 1) {
2228 mask.setActiveState(iter.getCoord() & ~1u,
true);
2234 for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
2235 for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
2236 for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
2237 if (!mask.isValueOn(ijk) &
isNonManifold(mDistAcc, ijk, mIsovalue, dim)) {
2238 mask.setActiveState(ijk,
true);
2245 gradients.setValuesOff();
2246 for (iter = mSignLeafs.leaf(n).cbeginValueOn(); iter; ++iter) {
2247 ijk = iter.getCoord();
2248 if(!mask.isValueOn(ijk & ~1u)) {
2251 gradients.setValueOn(iter.pos(), dir);
2257 for ( ; dim <= LeafDim; dim = dim << 1) {
2258 const unsigned coordMask = ~((dim << 1) - 1);
2259 for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
2260 for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
2261 for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
2263 adaptivity = adaptivityLeaf.getValue(ijk);
2265 if (mask.isValueOn(ijk) ||
isNonManifold(mDistAcc, ijk, mIsovalue, dim)
2266 || !
isMergable(gradients, ijk, dim, adaptivity)) {
2267 mask.setActiveState(ijk & coordMask,
true);
2289 mPolygonPool = &quadPool;
2297 mPolygonPool->
quad(mIdx) = verts;
2327 mPolygonPool = &polygonPool;
2337 if (verts[0] != verts[1] && verts[0] != verts[2] && verts[0] != verts[3]
2338 && verts[1] != verts[2] && verts[1] != verts[3] && verts[2] != verts[3]) {
2339 mPolygonPool->
quadFlags(mQuadIdx) = flags;
2340 addQuad(verts, reverse);
2342 verts[0] == verts[3] &&
2343 verts[1] != verts[2] &&
2344 verts[1] != verts[0] &&
2345 verts[2] != verts[0]) {
2347 addTriangle(verts[0], verts[1], verts[2], reverse);
2349 verts[1] == verts[2] &&
2350 verts[0] != verts[3] &&
2351 verts[0] != verts[1] &&
2352 verts[3] != verts[1]) {
2354 addTriangle(verts[0], verts[1], verts[3], reverse);
2356 verts[0] == verts[1] &&
2357 verts[2] != verts[3] &&
2358 verts[2] != verts[0] &&
2359 verts[3] != verts[0]) {
2361 addTriangle(verts[0], verts[2], verts[3], reverse);
2363 verts[2] == verts[3] &&
2364 verts[0] != verts[1] &&
2365 verts[0] != verts[2] &&
2366 verts[1] != verts[2]) {
2368 addTriangle(verts[0], verts[1], verts[2], reverse);
2375 mPolygonPool->
trimQuads(mQuadIdx,
true);
2381 void addQuad(
const Vec4I& verts,
bool reverse)
2384 mPolygonPool->
quad(mQuadIdx) = verts;
2386 Vec4I& quad = mPolygonPool->
quad(mQuadIdx);
2395 void addTriangle(
unsigned v0,
unsigned v1,
unsigned v2,
bool reverse)
2411 size_t mQuadIdx, mTriangleIdx;
2412 PolygonPool *mPolygonPool;
2416 template<
typename SignAccT,
typename IdxAccT,
typename PrimBuilder>
2419 const SignAccT& signAcc,
const IdxAccT& idxAcc, PrimBuilder& mesher,
Index32 pointListSize)
2421 const Index32 v0 = idxAcc.getValue(ijk);
2428 const bool isInside = flags &
INSIDE;
2435 if (flags &
XEDGE) {
2437 quad[0] = v0 + offsets[0];
2441 coord[1] = ijk[1] - 1;
2444 quad[1] = idxAcc.getValue(coord);
2445 cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2446 if (sEdgeGroupTable[cell][0] > 1) {
2447 tmpIdx = quad[1] +
Index32(sEdgeGroupTable[cell][5] - 1);
2448 if (tmpIdx < pointListSize) quad[1] = tmpIdx;
2454 quad[2] = idxAcc.getValue(coord);
2455 cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2456 if (sEdgeGroupTable[cell][0] > 1) {
2457 tmpIdx = quad[2] +
Index32(sEdgeGroupTable[cell][7] - 1);
2458 if (tmpIdx < pointListSize) quad[2] = tmpIdx;
2464 quad[3] = idxAcc.getValue(coord);
2465 cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2466 if (sEdgeGroupTable[cell][0] > 1) {
2467 tmpIdx = quad[3] +
Index32(sEdgeGroupTable[cell][3] - 1);
2468 if (tmpIdx < pointListSize) quad[3] = tmpIdx;
2473 mesher.addPrim(quad, isInside, tag[
bool(refFlags & XEDGE)]);
2478 if (flags &
YEDGE) {
2480 quad[0] = v0 + offsets[1];
2485 coord[2] = ijk[2] - 1;
2487 quad[1] = idxAcc.getValue(coord);
2488 cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2489 if (sEdgeGroupTable[cell][0] > 1) {
2490 tmpIdx = quad[1] +
Index32(sEdgeGroupTable[cell][12] - 1);
2491 if (tmpIdx < pointListSize) quad[1] = tmpIdx;
2497 quad[2] = idxAcc.getValue(coord);
2498 cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2499 if (sEdgeGroupTable[cell][0] > 1) {
2500 tmpIdx = quad[2] +
Index32(sEdgeGroupTable[cell][11] - 1);
2501 if (tmpIdx < pointListSize) quad[2] = tmpIdx;
2507 quad[3] = idxAcc.getValue(coord);
2508 cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2509 if (sEdgeGroupTable[cell][0] > 1) {
2510 tmpIdx = quad[3] +
Index32(sEdgeGroupTable[cell][10] - 1);
2511 if (tmpIdx < pointListSize) quad[3] = tmpIdx;
2516 mesher.addPrim(quad, isInside, tag[
bool(refFlags & YEDGE)]);
2520 if (flags &
ZEDGE) {
2522 quad[0] = v0 + offsets[2];
2526 coord[1] = ijk[1] - 1;
2529 quad[1] = idxAcc.getValue(coord);
2530 cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2531 if (sEdgeGroupTable[cell][0] > 1) {
2532 tmpIdx = quad[1] +
Index32(sEdgeGroupTable[cell][8] - 1);
2533 if (tmpIdx < pointListSize) quad[1] = tmpIdx;
2539 quad[2] = idxAcc.getValue(coord);
2540 cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2541 if (sEdgeGroupTable[cell][0] > 1) {
2542 tmpIdx = quad[2] +
Index32(sEdgeGroupTable[cell][6] - 1);
2543 if (tmpIdx < pointListSize) quad[2] = tmpIdx;
2549 quad[3] = idxAcc.getValue(coord);
2550 cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2551 if (sEdgeGroupTable[cell][0] > 1) {
2552 tmpIdx = quad[3] +
Index32(sEdgeGroupTable[cell][2] - 1);
2553 if (tmpIdx < pointListSize) quad[3] = tmpIdx;
2558 mesher.addPrim(quad, !isInside, tag[
bool(refFlags & ZEDGE)]);
2567 template<
typename LeafManagerT,
typename PrimBuilder>
2571 typedef typename LeafManagerT::TreeType::template ValueConverter<int>::Type
IntTreeT;
2572 typedef typename LeafManagerT::TreeType::template ValueConverter<Int16>::Type
Int16TreeT;
2580 GenPolygons(
const LeafManagerT& signLeafs,
const Int16TreeT& signTree,
2581 const IntTreeT& idxTree, PolygonPoolList& polygons,
Index32 pointListSize);
2583 void run(
bool threaded =
true);
2591 void operator()(
const tbb::blocked_range<size_t>&)
const;
2594 const LeafManagerT& mSignLeafs;
2595 const Int16TreeT& mSignTree;
2596 const IntTreeT& mIdxTree;
2597 const PolygonPoolList& mPolygonPoolList;
2600 const Int16TreeT *mRefSignTree;
2604 template<
typename LeafManagerT,
typename PrimBuilder>
2608 : mSignLeafs(signLeafs)
2609 , mSignTree(signTree)
2611 , mPolygonPoolList(polygons)
2612 , mPointListSize(pointListSize)
2613 , mRefSignTree(NULL)
2617 template<
typename LeafManagerT,
typename PrimBuilder>
2621 if (threaded) tbb::parallel_for(mSignLeafs.getRange(), *
this);
2622 else (*
this)(mSignLeafs.getRange());
2625 template<
typename LeafManagerT,
typename PrimBuilder>
2628 const tbb::blocked_range<size_t>& range)
const
2630 typename LeafManagerT::TreeType::LeafNodeType::ValueOnCIter iter;
2641 boost::scoped_ptr<Int16AccessorT> refSignAcc;
2642 if (mRefSignTree) refSignAcc.reset(
new Int16AccessorT(*mRefSignTree));
2645 for (
size_t n = range.begin(); n != range.end(); ++n) {
2647 origin = mSignLeafs.leaf(n).origin();
2651 iter = mSignLeafs.leaf(n).cbeginValueOn();
2652 for (; iter; ++iter) {
2653 if (iter.getValue() &
XEDGE) ++edgeCount;
2654 if (iter.getValue() &
YEDGE) ++edgeCount;
2655 if (iter.getValue() &
ZEDGE) ++edgeCount;
2658 if(edgeCount == 0)
continue;
2660 mesher.init(edgeCount, mPolygonPoolList[n]);
2662 const typename Int16TreeT::LeafNodeType *signleafPt = signAcc.
probeConstLeaf(origin);
2663 const typename IntTreeT::LeafNodeType *idxLeafPt = idxAcc.
probeConstLeaf(origin);
2665 if (!signleafPt || !idxLeafPt)
continue;
2668 const typename Int16TreeT::LeafNodeType *refSignLeafPt = NULL;
2669 if (refSignAcc) refSignLeafPt = refSignAcc->probeConstLeaf(origin);
2673 iter = mSignLeafs.leaf(n).cbeginValueOn();
2674 for (; iter; ++iter) {
2675 ijk = iter.getCoord();
2677 Int16 flags = iter.getValue();
2679 if (!(flags & 0xE00))
continue;
2682 if (refSignLeafPt) {
2683 refFlags = refSignLeafPt->getValue(iter.pos());
2690 const unsigned char cell = uint8_t(
SIGNS & flags);
2692 if (sEdgeGroupTable[cell][0] > 1) {
2693 offsets[0] = (sEdgeGroupTable[cell][1] - 1);
2694 offsets[1] = (sEdgeGroupTable[cell][9] - 1);
2695 offsets[2] = (sEdgeGroupTable[cell][4] - 1);
2698 if (ijk[0] > origin[0] && ijk[1] > origin[1] && ijk[2] > origin[2]) {
2700 *signleafPt, *idxLeafPt, mesher, mPointListSize);
2703 signAcc, idxAcc, mesher, mPointListSize);
2719 PartOp(
size_t leafCount,
size_t partitions,
size_t activePart)
2721 size_t leafSegments = leafCount / partitions;
2722 mStart = leafSegments * activePart;
2723 mEnd = activePart >= (partitions - 1) ? leafCount : mStart + leafSegments;
2726 template <
typename LeafNodeType>
2729 if (leafIndex < mStart || leafIndex >= mEnd) leaf.setValuesOff();
2733 size_t mStart, mEnd;
2740 template<
typename SrcTreeT>
2745 typedef typename SrcTreeT::template ValueConverter<bool>::Type
BoolTreeT;
2751 PartGen(
const LeafManagerT& leafs,
size_t partitions,
size_t activePart);
2753 void run(
bool threaded =
true);
2755 BoolTreeT&
tree() {
return mTree; }
2761 void operator()(
const tbb::blocked_range<size_t>&);
2765 const LeafManagerT& mLeafManager;
2767 size_t mStart, mEnd;
2770 template<
typename SrcTreeT>
2772 : mLeafManager(leafs)
2778 size_t leafSegments = leafCount / partitions;
2779 mStart = leafSegments * activePart;
2780 mEnd = activePart >= (partitions - 1) ? leafCount : mStart + leafSegments;
2783 template<
typename SrcTreeT>
2785 : mLeafManager(rhs.mLeafManager)
2787 , mStart(rhs.mStart)
2793 template<
typename SrcTreeT>
2797 if (threaded) tbb::parallel_reduce(mLeafManager.getRange(), *
this);
2798 else (*
this)(mLeafManager.getRange());
2802 template<
typename SrcTreeT>
2809 typedef typename BoolTreeT::LeafNodeType BoolLeafT;
2810 typename SrcTreeT::LeafNodeType::ValueOnCIter iter;
2812 for (
size_t n = range.begin(); n != range.end(); ++n) {
2813 if (n < mStart || n >= mEnd)
continue;
2814 BoolLeafT* leaf = acc.
touchLeaf(mLeafManager.leaf(n).origin());
2815 leaf->topologyUnion(mLeafManager.leaf(n));
2823 template<
typename TreeT,
typename LeafManagerT>
2827 typedef typename TreeT::template ValueConverter<bool>::Type
BoolTreeT;
2831 GenSeamMask(
const LeafManagerT& leafs,
const TreeT& tree);
2833 void run(
bool threaded =
true);
2835 BoolTreeT&
mask() {
return mMaskTree; }
2840 void operator()(
const tbb::blocked_range<size_t>&);
2845 const LeafManagerT& mLeafManager;
2848 BoolTreeT mMaskTree;
2852 template<
typename TreeT,
typename LeafManagerT>
2854 : mLeafManager(leafs)
2861 template<
typename TreeT,
typename LeafManagerT>
2863 : mLeafManager(rhs.mLeafManager)
2870 template<
typename TreeT,
typename LeafManagerT>
2874 if (threaded) tbb::parallel_reduce(mLeafManager.getRange(), *
this);
2875 else (*
this)(mLeafManager.getRange());
2879 template<
typename TreeT,
typename LeafManagerT>
2887 typename LeafManagerT::TreeType::LeafNodeType::ValueOnCIter it;
2889 for (
size_t n = range.begin(); n != range.end(); ++n) {
2891 it = mLeafManager.leaf(n).cbeginValueOn();
2895 ijk = it.getCoord();
2897 unsigned char rhsSigns = uint8_t(acc.
getValue(ijk) &
SIGNS);
2899 if (sEdgeGroupTable[rhsSigns][0] > 0) {
2900 unsigned char lhsSigns = uint8_t(it.getValue() &
SIGNS);
2901 if (rhsSigns != lhsSigns) {
2913 template<
typename TreeT>
2921 template <
typename LeafNodeType>
2924 const typename TreeT::LeafNodeType *maskLeaf =
2927 if (!maskLeaf)
return;
2929 typename LeafNodeType::ValueOnIter it = leaf.beginValueOn();
2933 if (maskLeaf->isValueOn(it.pos())) {
2934 it.setValue(it.getValue() |
SEAM);
2945 template<
typename BoolTreeT>
2950 MaskEdges(
const BoolTreeT& valueMask) : mMaskAcc(valueMask) {}
2952 template <
typename LeafNodeType>
2955 typename LeafNodeType::ValueOnIter it = leaf.beginValueOn();
2957 const typename BoolTreeT::LeafNodeType * maskLeaf =
2962 if (!maskLeaf->isValueOn(it.pos())) {
2963 it.setValue(0x1FF & it.getValue());
2968 it.setValue(0x1FF & it.getValue());
2974 BoolAccessorT mMaskAcc;
2984 std::vector<unsigned char>& usedPointMask)
2985 : mPolygons(polygons)
2986 , mPolyListCount(polyListCount)
2987 , mUsedPointMask(usedPointMask)
2991 void run(
bool threaded =
true)
2994 tbb::parallel_for(tbb::blocked_range<size_t>(0, mPolyListCount), *
this);
2996 (*this)(tbb::blocked_range<size_t>(0, mPolyListCount));
3006 for (
size_t n = range.begin(); n != range.end(); ++n) {
3008 for (
size_t i = 0; i < polygons.
numQuads(); ++i) {
3010 mUsedPointMask[quad[0]] = 1;
3011 mUsedPointMask[quad[1]] = 1;
3012 mUsedPointMask[quad[2]] = 1;
3013 mUsedPointMask[quad[3]] = 1;
3018 mUsedPointMask[triangle[0]] = 1;
3019 mUsedPointMask[triangle[1]] = 1;
3020 mUsedPointMask[triangle[2]] = 1;
3027 const PolygonPoolList& mPolygons;
3028 size_t mPolyListCount;
3029 std::vector<unsigned char>& mUsedPointMask;
3038 size_t polyListCount,
const std::vector<unsigned>& indexMap)
3039 : mPolygons(polygons)
3040 , mPolyListCount(polyListCount)
3041 , mIndexMap(indexMap)
3045 void run(
bool threaded =
true)
3048 tbb::parallel_for(tbb::blocked_range<size_t>(0, mPolyListCount), *
this);
3050 (*this)(tbb::blocked_range<size_t>(0, mPolyListCount));
3058 for (
size_t n = range.begin(); n != range.end(); ++n) {
3060 for (
size_t i = 0; i < polygons.
numQuads(); ++i) {
3062 quad[0] = mIndexMap[quad[0]];
3063 quad[1] = mIndexMap[quad[1]];
3064 quad[2] = mIndexMap[quad[2]];
3065 quad[3] = mIndexMap[quad[3]];
3070 triangle[0] = mIndexMap[triangle[0]];
3071 triangle[1] = mIndexMap[triangle[1]];
3072 triangle[2] = mIndexMap[triangle[2]];
3079 PolygonPoolList& mPolygons;
3080 size_t mPolyListCount;
3081 const std::vector<unsigned>& mIndexMap;
3092 const PointList& oldPointList,
3093 const std::vector<unsigned>& indexMap,
3094 const std::vector<unsigned char>& usedPointMask)
3095 : mNewPointList(newPointList)
3096 , mOldPointList(oldPointList)
3097 , mIndexMap(indexMap)
3098 , mUsedPointMask(usedPointMask)
3102 void run(
bool threaded =
true)
3105 tbb::parallel_for(tbb::blocked_range<size_t>(0, mIndexMap.size()), *
this);
3107 (*this)(tbb::blocked_range<size_t>(0, mIndexMap.size()));
3115 for (
size_t n = range.begin(); n != range.end(); ++n) {
3116 if (mUsedPointMask[n]) {
3117 const size_t index = mIndexMap[n];
3118 mNewPointList.get()[index] = mOldPointList[n];
3125 const PointList& mOldPointList;
3126 const std::vector<unsigned>& mIndexMap;
3127 const std::vector<unsigned char>& mUsedPointMask;
3134 template<
typename SrcTreeT>
3139 typedef typename SrcTreeT::template ValueConverter<bool>::Type
BoolTreeT;
3151 void run(
bool threaded =
true);
3153 BoolTreeT&
tree() {
return mTree; }
3160 void operator()(
const tbb::blocked_range<size_t>&);
3166 const BoolGridT& mMask;
3167 const LeafManagerT& mLeafManager;
3174 template<
typename SrcTreeT>
3178 , mLeafManager(srcLeafs)
3179 , mSrcXForm(srcXForm)
3180 , mInvertMask(invertMask)
3186 template<
typename SrcTreeT>
3189 , mLeafManager(rhs.mLeafManager)
3190 , mSrcXForm(rhs.mSrcXForm)
3191 , mInvertMask(rhs.mInvertMask)
3197 template<
typename SrcTreeT>
3202 tbb::parallel_reduce(mLeafManager.getRange(), *
this);
3204 (*this)(mLeafManager.getRange());
3209 template<
typename SrcTreeT>
3215 typedef typename BoolTreeT::LeafNodeType BoolLeafT;
3220 typename SrcTreeT::LeafNodeType::ValueOnCIter iter;
3221 for (
size_t n = range.begin(); n != range.end(); ++n) {
3223 ijk = mLeafManager.leaf(n).origin();
3224 BoolLeafT* leaf =
new BoolLeafT(ijk,
false);
3225 bool addLeaf =
false;
3227 if (maskXForm == mSrcXForm) {
3229 const BoolLeafT* maskLeaf = maskAcc.probeConstLeaf(ijk);
3233 for (iter = mLeafManager.leaf(n).cbeginValueOn(); iter; ++iter) {
3234 Index pos = iter.pos();
3235 if(maskLeaf->isValueOn(pos) != mInvertMask) {
3236 leaf->setValueOn(pos);
3241 }
else if (maskAcc.isValueOn(ijk) != mInvertMask) {
3242 leaf->topologyUnion(mLeafManager.leaf(n));
3247 for (iter = mLeafManager.leaf(n).cbeginValueOn(); iter; ++iter) {
3248 ijk = iter.getCoord();
3249 xyz = maskXForm.
worldToIndex(mSrcXForm.indexToWorld(ijk));
3251 leaf->setValueOn(iter.pos());
3257 if (addLeaf) acc.
addLeaf(leaf);
3266 template<
typename SrcTreeT>
3270 typedef typename SrcTreeT::template ValueConverter<int>::Type
IntTreeT;
3271 typedef typename SrcTreeT::template ValueConverter<bool>::Type
BoolTreeT;
3276 GenBoundaryMask(
const LeafManagerT& leafs,
const BoolTreeT&,
const IntTreeT&);
3278 void run(
bool threaded =
true);
3280 BoolTreeT&
tree() {
return mTree; }
3285 void operator()(
const tbb::blocked_range<size_t>&);
3292 bool neighboringLeaf(
const Coord&,
const IntTreeAccessorT&)
const;
3294 const LeafManagerT& mLeafManager;
3295 const BoolTreeT& mMaskTree;
3296 const IntTreeT& mIdxTree;
3298 CoordBBox mLeafBBox;
3302 template<
typename SrcTreeT>
3305 : mLeafManager(leafs)
3306 , mMaskTree(maskTree)
3310 mIdxTree.evalLeafBoundingBox(mLeafBBox);
3311 mLeafBBox.expand(IntTreeT::LeafNodeType::DIM);
3315 template<
typename SrcTreeT>
3317 : mLeafManager(rhs.mLeafManager)
3318 , mMaskTree(rhs.mMaskTree)
3319 , mIdxTree(rhs.mIdxTree)
3321 , mLeafBBox(rhs.mLeafBBox)
3326 template<
typename SrcTreeT>
3331 tbb::parallel_reduce(mLeafManager.getRange(), *
this);
3333 (*this)(mLeafManager.getRange());
3338 template<
typename SrcTreeT>
3342 if (acc.probeConstLeaf(ijk))
return true;
3344 const int dim = IntTreeT::LeafNodeType::DIM;
3347 if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1], ijk[2])))
return true;
3348 if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1], ijk[2])))
return true;
3349 if (acc.probeConstLeaf(Coord(ijk[0], ijk[1] + dim, ijk[2])))
return true;
3350 if (acc.probeConstLeaf(Coord(ijk[0], ijk[1] - dim, ijk[2])))
return true;
3351 if (acc.probeConstLeaf(Coord(ijk[0], ijk[1], ijk[2] + dim)))
return true;
3352 if (acc.probeConstLeaf(Coord(ijk[0], ijk[1], ijk[2] - dim)))
return true;
3355 if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1], ijk[2] - dim)))
return true;
3356 if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1], ijk[2] - dim)))
return true;
3357 if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1], ijk[2] + dim)))
return true;
3358 if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1], ijk[2] + dim)))
return true;
3359 if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1] + dim, ijk[2])))
return true;
3360 if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1] + dim, ijk[2])))
return true;
3361 if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1] - dim, ijk[2])))
return true;
3362 if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1] - dim, ijk[2])))
return true;
3363 if (acc.probeConstLeaf(Coord(ijk[0], ijk[1] - dim, ijk[2] + dim)))
return true;
3364 if (acc.probeConstLeaf(Coord(ijk[0], ijk[1] - dim, ijk[2] - dim)))
return true;
3365 if (acc.probeConstLeaf(Coord(ijk[0], ijk[1] + dim, ijk[2] + dim)))
return true;
3366 if (acc.probeConstLeaf(Coord(ijk[0], ijk[1] + dim, ijk[2] - dim)))
return true;
3369 if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1] - dim, ijk[2] - dim)))
return true;
3370 if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1] - dim, ijk[2] + dim)))
return true;
3371 if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1] - dim, ijk[2] + dim)))
return true;
3372 if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1] - dim, ijk[2] - dim)))
return true;
3373 if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1] + dim, ijk[2] - dim)))
return true;
3374 if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1] + dim, ijk[2] + dim)))
return true;
3375 if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1] + dim, ijk[2] + dim)))
return true;
3376 if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1] + dim, ijk[2] - dim)))
return true;
3382 template<
typename SrcTreeT>
3391 typename SrcTreeT::LeafNodeType::ValueOnCIter iter;
3393 for (
size_t n = range.begin(); n != range.end(); ++n) {
3395 const typename SrcTreeT::LeafNodeType&
3396 leaf = mLeafManager.leaf(n);
3398 ijk = leaf.origin();
3400 if (!mLeafBBox.isInside(ijk) || !neighboringLeaf(ijk, idxAcc))
continue;
3402 const typename BoolTreeT::LeafNodeType*
3405 if (!maskLeaf || !leaf.hasSameTopology(maskLeaf)) {
3406 acc.
touchLeaf(ijk)->topologyUnion(leaf);
3415 template<
typename TreeT>
3419 typedef typename TreeT::template ValueConverter<bool>::Type
BoolTreeT;
3425 GenTileMask(
const std::vector<Vec4i>& tiles,
const TreeT& distTree, ValueT iso);
3427 void run(
bool threaded =
true);
3429 BoolTreeT&
tree() {
return mTree; }
3434 void operator()(
const tbb::blocked_range<size_t>&);
3439 const std::vector<Vec4i>& mTiles;
3440 const TreeT& mDistTree;
3447 template<
typename TreeT>
3449 const std::vector<Vec4i>& tiles,
const TreeT& distTree,
ValueT iso)
3451 , mDistTree(distTree)
3458 template<
typename TreeT>
3460 : mTiles(rhs.mTiles)
3461 , mDistTree(rhs.mDistTree)
3462 , mIsovalue(rhs.mIsovalue)
3468 template<
typename TreeT>
3472 if (threaded) tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mTiles.size()), *
this);
3473 else (*
this)(tbb::blocked_range<size_t>(0, mTiles.size()));
3477 template<
typename TreeT>
3482 CoordBBox region, bbox;
3484 bool processRegion =
true;
3488 for (
size_t n = range.begin(); n != range.end(); ++n) {
3490 const Vec4i& tile = mTiles[n];
3492 bbox.min()[0] = tile[0];
3493 bbox.min()[1] = tile[1];
3494 bbox.min()[2] = tile[2];
3496 bbox.max() = bbox.min();
3497 bbox.max().offset(tile[3]);
3499 const bool thisInside = (distAcc.
getValue(bbox.min()) < mIsovalue);
3508 processRegion =
true;
3510 processRegion = thisInside != (distAcc.
getValue(nijk) < mIsovalue);
3514 if (processRegion) {
3516 region.min()[0] = region.max()[0] = ijk[0];
3517 mTree.fill(region,
true);
3524 processRegion =
true;
3526 processRegion = !distAcc.
probeValue(ijk, value) && thisInside != (value < mIsovalue);
3529 if (processRegion) {
3531 region.min()[0] = region.max()[0] = ijk[0];
3532 mTree.fill(region,
true);
3542 processRegion =
true;
3544 processRegion = thisInside != (distAcc.
getValue(nijk) < mIsovalue);
3547 if (processRegion) {
3549 region.min()[1] = region.max()[1] = ijk[1];
3550 mTree.fill(region,
true);
3557 processRegion =
true;
3559 processRegion = !distAcc.
probeValue(ijk, value) && thisInside != (value < mIsovalue);
3562 if (processRegion) {
3564 region.min()[1] = region.max()[1] = ijk[1];
3565 mTree.fill(region,
true);
3575 processRegion =
true;
3577 processRegion = thisInside != (distAcc.
getValue(nijk) < mIsovalue);
3580 if (processRegion) {
3582 region.min()[2] = region.max()[2] = ijk[2];
3583 mTree.fill(region,
true);
3589 processRegion =
true;
3591 processRegion = !distAcc.
probeValue(ijk, value) && thisInside != (value < mIsovalue);
3594 if (processRegion) {
3596 region.min()[2] = region.max()[2] = ijk[2];
3597 mTree.fill(region,
true);
3605 processRegion =
true;
3607 processRegion = !distAcc.
probeValue(ijk, value) && thisInside != (value < mIsovalue);
3610 if (processRegion) {
3612 region.min()[1] = region.max()[1] = ijk[1];
3613 region.min()[2] = region.max()[2] = ijk[2];
3614 mTree.fill(region,
true);
3622 processRegion =
true;
3624 processRegion = !distAcc.
probeValue(ijk, value) && thisInside != (value < mIsovalue);
3627 if (processRegion) {
3629 region.min()[1] = region.max()[1] = ijk[1];
3630 region.min()[0] = region.max()[0] = ijk[0];
3631 mTree.fill(region,
true);
3638 processRegion =
true;
3640 processRegion = !distAcc.
probeValue(ijk, value) && thisInside != (value < mIsovalue);
3643 if (processRegion) {
3645 region.min()[2] = region.max()[2] = ijk[2];
3646 region.min()[0] = region.max()[0] = ijk[0];
3647 mTree.fill(region,
true);
3656 template<
class DistTreeT,
class SignTreeT,
class IdxTreeT>
3658 tileData(
const DistTreeT& distTree, SignTreeT& signTree, IdxTreeT& idxTree,
double iso)
3660 typename DistTreeT::ValueOnCIter tileIter(distTree);
3661 tileIter.setMaxDepth(DistTreeT::ValueOnCIter::LEAF_DEPTH - 1);
3663 if (!tileIter)
return;
3665 size_t tileCount = 0;
3666 for ( ; tileIter; ++tileIter) {
3670 std::vector<Vec4i> tiles(tileCount);
3673 tileIter = distTree.cbeginValueOn();
3674 tileIter.setMaxDepth(DistTreeT::ValueOnCIter::LEAF_DEPTH - 1);
3677 for (; tileIter; ++tileIter) {
3678 Vec4i& tile = tiles[tileCount++];
3679 tileIter.getBoundingBox(bbox);
3680 tile[0] = bbox.min()[0];
3681 tile[1] = bbox.min()[1];
3682 tile[2] = bbox.min()[2];
3683 tile[3] = bbox.max()[0] - bbox.min()[0];
3686 typename DistTreeT::ValueType isovalue =
typename DistTreeT::ValueType(iso);
3691 typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
3694 BoolLeafManagerT leafs(tileMask.
tree());
3700 signTree.merge(*op.signTree());
3701 idxTree.merge(*op.idxTree());
3713 : mPointsIn(pointsIn) , mPointsOut(pointsOut)
3719 for (
size_t n = range.begin(); n < range.end(); ++n) {
3720 mPointsOut[n] = mPointsIn[n];
3725 const PointList& mPointsIn;
3726 std::vector<Vec3s>& mPointsOut;
3731 template <
typename LeafManagerT>
3735 double interiorWidth = 0.0, exteriorWidth = 0.0;
3737 typename LeafManagerT::TreeType::LeafNodeType::ValueOffCIter it;
3738 bool foundInterior =
false, foundExterior =
false;
3739 for (
size_t n = 0, N = leafs.leafCount(); n < N; ++n) {
3741 for (it = leafs.leaf(n).cbeginValueOff(); it; ++it) {
3742 double value = double(it.getValue());
3744 interiorWidth = value;
3745 foundInterior =
true;
3746 }
else if (value > 0.0) {
3747 exteriorWidth = value;
3748 foundExterior =
true;
3751 if (foundInterior && foundExterior)
break;
3754 if (foundInterior && foundExterior)
break;
3759 double minDist =
std::min(std::abs(interiorWidth - iso), std::abs(exteriorWidth - iso));
3760 return !(minDist > (2.0 * voxelSize));
3777 , mTriangleFlags(NULL)
3784 : mNumQuads(numQuads)
3785 , mNumTriangles(numTriangles)
3788 , mQuadFlags(new char[mNumQuads])
3789 , mTriangleFlags(new char[mNumTriangles])
3800 for (
size_t i = 0; i < mNumQuads; ++i) {
3801 mQuads[i] = rhs.mQuads[i];
3802 mQuadFlags[i] = rhs.mQuadFlags[i];
3805 for (
size_t i = 0; i < mNumTriangles; ++i) {
3806 mTriangles[i] = rhs.mTriangles[i];
3807 mTriangleFlags[i] = rhs.mTriangleFlags[i];
3817 mQuadFlags.reset(
new char[mNumQuads]);
3826 mQuadFlags.reset(NULL);
3833 mNumTriangles = size;
3835 mTriangleFlags.reset(
new char[mNumTriangles]);
3843 mTriangles.reset(NULL);
3844 mTriangleFlags.reset(NULL);
3851 if (!(n < mNumQuads))
return false;
3859 boost::scoped_array<openvdb::Vec4I> quads(
new openvdb::Vec4I[n]);
3860 boost::scoped_array<char> flags(
new char[n]);
3862 for (
size_t i = 0; i < n; ++i) {
3863 quads[i] = mQuads[i];
3864 flags[i] = mQuadFlags[i];
3868 mQuadFlags.swap(flags);
3880 if (!(n < mNumTriangles))
return false;
3885 mTriangles.reset(NULL);
3888 boost::scoped_array<openvdb::Vec3I> triangles(
new openvdb::Vec3I[n]);
3889 boost::scoped_array<char> flags(
new char[n]);
3891 for (
size_t i = 0; i < n; ++i) {
3892 triangles[i] = mTriangles[i];
3893 flags[i] = mTriangleFlags[i];
3896 mTriangles.swap(triangles);
3897 mTriangleFlags.swap(flags);
3913 , mSeamPointListSize(0)
3914 , mPolygonPoolListSize(0)
3915 , mIsovalue(isovalue)
3916 , mPrimAdaptivity(adaptivity)
3917 , mSecAdaptivity(0.0)
3919 , mSurfaceMaskGrid(
GridBase::ConstPtr())
3920 , mAdaptivityGrid(
GridBase::ConstPtr())
3921 , mAdaptivityMaskTree(
TreeBase::ConstPtr())
3924 , mInvertSurfaceMask(false)
3927 , mQuantizedSeamPoints(NULL)
3940 inline const size_t&
3943 return mPointListSize;
3947 inline PolygonPoolList&
3954 inline const PolygonPoolList&
3961 inline const size_t&
3964 return mPolygonPoolListSize;
3972 mSecAdaptivity = secAdaptivity;
3977 mSeamPointListSize = 0;
3978 mQuantizedSeamPoints.reset(NULL);
3985 mSurfaceMaskGrid = mask;
3986 mInvertSurfaceMask = invertMask;
3993 mAdaptivityGrid = grid;
4000 mAdaptivityMaskTree = tree;
4007 mPartitions =
std::max(partitions,
unsigned(1));
4008 mActivePart =
std::min(activePart, mPartitions-1);
4012 inline std::vector<unsigned char>&
4019 inline const std::vector<unsigned char>&
4026 template<
typename Gr
idT>
4030 typedef typename GridT::TreeType DistTreeT;
4032 typedef typename DistTreeT::ValueType DistValueT;
4034 typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
4038 typedef typename DistTreeT::template ValueConverter<Int16>::Type Int16TreeT;
4041 typedef typename DistTreeT::template ValueConverter<int>::Type IntTreeT;
4042 typedef typename DistTreeT::template ValueConverter<float>::Type FloatTreeT;
4046 const openvdb::math::Transform& transform = distGrid.
transform();
4047 const DistTreeT& distTree = distGrid.tree();
4048 const DistValueT isovalue = DistValueT(mIsovalue);
4050 typename Int16TreeT::Ptr signTreePt;
4051 typename IntTreeT::Ptr idxTreePt;
4052 typename BoolTreeT::Ptr pointMask;
4054 BoolTreeT valueMask(
false), seamMask(
false);
4055 const bool adaptive = mPrimAdaptivity > 1e-7 || mSecAdaptivity > 1e-7;
4056 bool maskEdges =
false;
4059 const BoolGridT * surfaceMask = NULL;
4060 if (mSurfaceMaskGrid && mSurfaceMaskGrid->type() == BoolGridT::gridType()) {
4061 surfaceMask =
static_cast<const BoolGridT*
>(mSurfaceMaskGrid.get());
4064 const FloatGridT * adaptivityField = NULL;
4065 if (mAdaptivityGrid && mAdaptivityGrid->type() == FloatGridT::gridType()) {
4066 adaptivityField =
static_cast<const FloatGridT*
>(mAdaptivityGrid.get());
4069 if (mAdaptivityMaskTree && mAdaptivityMaskTree->type() == BoolTreeT::treeType()) {
4070 const BoolTreeT *adaptivityMaskPt =
4071 static_cast<const BoolTreeT*
>(mAdaptivityMaskTree.get());
4072 seamMask.topologyUnion(*adaptivityMaskPt);
4078 DistLeafManagerT distLeafs(distTree);
4081 bool padActiveVoxels =
false;
4085 padActiveVoxels =
true;
4088 mIsovalue, transform.voxelSize()[0]);
4092 if (!padActiveVoxels) {
4094 distTree.evalActiveVoxelDim(dim);
4096 if (maxDim < 1000) {
4097 padActiveVoxels =
true;
4102 if (surfaceMask || mPartitions > 1) {
4110 *surfaceMask, distLeafs, transform, mInvertSurfaceMask);
4112 valueMask.merge(masking.
tree());
4115 if (mPartitions > 1) {
4125 valueMask.merge(partitioner.
tree());
4130 BoolLeafManagerT leafs(valueMask);
4133 signDataOp(distTree, leafs, isovalue);
4136 signTreePt = signDataOp.
signTree();
4137 idxTreePt = signDataOp.
idxTree();
4144 BoolLeafManagerT bleafs(boundary.
tree());
4147 signDataOp(distTree, bleafs, isovalue);
4150 signTreePt->merge(*signDataOp.signTree());
4151 idxTreePt->merge(*signDataOp.idxTree());
4157 if (padActiveVoxels) {
4159 BoolTreeT regionMask(
false);
4160 regionMask.topologyUnion(distTree);
4163 BoolLeafManagerT leafs(regionMask);
4166 signDataOp(distTree, leafs, isovalue);
4169 signTreePt = signDataOp.
signTree();
4170 idxTreePt = signDataOp.
idxTree();
4174 signDataOp(distTree, distLeafs, isovalue);
4177 signTreePt = signDataOp.
signTree();
4178 idxTreePt = signDataOp.
idxTree();
4186 internal::tileData(distTree, *signTreePt, *idxTreePt, static_cast<double>(isovalue));
4189 Int16TreeT *refSignTreePt = NULL;
4190 IntTreeT *refIdxTreePt = NULL;
4191 const DistTreeT *refDistTreePt = NULL;
4193 if (mRefGrid && mRefGrid->type() == GridT::gridType()) {
4195 const GridT* refGrid =
static_cast<const GridT*
>(mRefGrid.get());
4196 refDistTreePt = &refGrid->tree();
4199 if (!mRefSignTree && !mRefIdxTree) {
4201 DistLeafManagerT refDistLeafs(*refDistTreePt);
4203 signDataOp(*refDistTreePt, refDistLeafs, isovalue);
4207 mRefSignTree = signDataOp.
signTree();
4208 mRefIdxTree = signDataOp.
idxTree();
4212 if (mRefSignTree && mRefIdxTree) {
4213 refSignTreePt =
static_cast<Int16TreeT*
>(mRefSignTree.get());
4214 refIdxTreePt =
static_cast<IntTreeT*
>(mRefIdxTree.get());
4220 Int16LeafManagerT signLeafs(*signTreePt);
4229 if (refSignTreePt) {
4236 seamMask.merge(seamOp.
mask());
4240 std::vector<size_t> regions(signLeafs.leafCount(), 0);
4241 if (regions.empty())
return;
4246 signLeafs, *signTreePt, distTree, *idxTreePt, isovalue, DistValueT(mPrimAdaptivity));
4248 if (adaptivityField) {
4252 if (refSignTreePt || mAdaptivityMaskTree) {
4256 if (refSignTreePt) {
4257 merge.
setRefData(refSignTreePt, DistValueT(mSecAdaptivity));
4272 for (
size_t n = 0, N = regions.size(); n < N; ++n) {
4274 regions[n] = mPointListSize;
4275 mPointListSize += tmp;
4282 mPointFlags.clear();
4285 if (refSignTreePt && refIdxTreePt) {
4287 if (mSeamPointListSize == 0) {
4289 std::vector<size_t> pointMap;
4292 Int16LeafManagerT refSignLeafs(*refSignTreePt);
4293 pointMap.resize(refSignLeafs.leafCount(), 0);
4298 for (
size_t n = 0, N = pointMap.size(); n < N; ++n) {
4300 pointMap[n] = mSeamPointListSize;
4301 mSeamPointListSize += tmp;
4305 if (!pointMap.empty() && mSeamPointListSize != 0) {
4307 mQuantizedSeamPoints.reset(
new uint32_t[mSeamPointListSize]);
4308 memset(mQuantizedSeamPoints.get(), 0,
sizeof(uint32_t) * mSeamPointListSize);
4312 IntLeafManagerT refIdxLeafs(*refIdxTreePt);
4317 if (mSeamPointListSize != 0) {
4319 distTree, *refSignTreePt, *refIdxTreePt, mQuantizedSeamPoints, mIsovalue));
4325 pointOp(signLeafs, distTree, *idxTreePt, mPoints, regions, transform, mIsovalue);
4328 if (mSeamPointListSize != 0) {
4329 mPointFlags.resize(mPointListSize);
4330 pointOp.
setRefData(refSignTreePt, refDistTreePt, refIdxTreePt,
4331 &mQuantizedSeamPoints, &mPointFlags);
4337 mPolygonPoolListSize = signLeafs.leafCount();
4338 mPolygons.reset(
new PolygonPool[mPolygonPoolListSize]);
4344 mesher(signLeafs, *signTreePt, *idxTreePt, mPolygons,
Index32(mPointListSize));
4352 mesher(signLeafs, *signTreePt, *idxTreePt, mPolygons,
Index32(mPointListSize));
4360 if ((surfaceMask || mPartitions > 1) && mPointListSize > 0) {
4363 std::vector<unsigned char> usedPointMask(mPointListSize, 0);
4369 std::vector<unsigned> indexMap(mPointListSize);
4370 size_t usedPointCount = 0;
4371 for (
size_t p = 0; p < mPointListSize; ++p) {
4372 if (usedPointMask[p]) indexMap[p] =
static_cast<unsigned>(usedPointCount++);
4375 if (usedPointCount < mPointListSize) {
4384 mPointListSize = usedPointCount;
4385 mPoints.reset(newPointList.release());
4396 if (refSignTreePt || refIdxTreePt || refDistTreePt) {
4397 std::vector<Vec3s> newPoints;
4399 for (
size_t n = 0; n < mPolygonPoolListSize; ++n) {
4403 std::vector<size_t> nonPlanarQuads;
4404 nonPlanarQuads.reserve(polygons.
numQuads());
4406 for (
size_t i = 0; i < polygons.
numQuads(); ++i) {
4414 const bool edgePoly = mPointFlags[quad[0]] || mPointFlags[quad[1]]
4415 || mPointFlags[quad[2]] || mPointFlags[quad[3]];
4417 if (!edgePoly)
continue;
4419 const Vec3s& p0 = mPoints[quad[0]];
4420 const Vec3s& p1 = mPoints[quad[1]];
4421 const Vec3s& p2 = mPoints[quad[2]];
4422 const Vec3s& p3 = mPoints[quad[3]];
4425 nonPlanarQuads.push_back(i);
4431 if (!nonPlanarQuads.empty()) {
4438 size_t triangleIdx = 0;
4439 for (
size_t i = 0; i < nonPlanarQuads.size(); ++i) {
4441 size_t& quadIdx = nonPlanarQuads[i];
4444 char& quadFlags = polygons.
quadFlags(quadIdx);
4447 Vec3s centroid = (mPoints[quad[0]] + mPoints[quad[1]] +
4448 mPoints[quad[2]] + mPoints[quad[3]]) * 0.25;
4450 size_t pointIdx = newPoints.size() + mPointListSize;
4452 newPoints.push_back(centroid);
4458 triangle[0] = quad[0];
4459 triangle[1] =
static_cast<unsigned>(pointIdx);
4460 triangle[2] = quad[3];
4464 if (mPointFlags[triangle[0]] || mPointFlags[triangle[2]]) {
4474 triangle[0] = quad[0];
4475 triangle[1] = quad[1];
4476 triangle[2] =
static_cast<unsigned>(pointIdx);
4480 if (mPointFlags[triangle[0]] || mPointFlags[triangle[1]]) {
4490 triangle[0] = quad[1];
4491 triangle[1] = quad[2];
4492 triangle[2] =
static_cast<unsigned>(pointIdx);
4496 if (mPointFlags[triangle[0]] || mPointFlags[triangle[1]]) {
4507 triangle[0] = quad[2];
4508 triangle[1] = quad[3];
4509 triangle[2] =
static_cast<unsigned>(pointIdx);
4513 if (mPointFlags[triangle[0]] || mPointFlags[triangle[1]]) {
4532 for (
size_t i = 0; i < polygons.
numQuads(); ++i) {
4536 tmpPolygons.
quad(quadIdx) = quad;
4543 polygons.
copy(tmpPolygons);
4549 if (!newPoints.empty()) {
4551 size_t newPointCount = newPoints.size() + mPointListSize;
4556 for (
size_t i = 0; i < mPointListSize; ++i) {
4557 newPointList.get()[i] = mPoints[i];
4560 for (
size_t i = mPointListSize; i < newPointCount; ++i) {
4561 newPointList.get()[i] = newPoints[i - mPointListSize];
4564 mPointListSize = newPointCount;
4565 mPoints.reset(newPointList.release());
4566 mPointFlags.resize(mPointListSize, 0);
4576 template<
typename Gr
idType>
4577 inline typename boost::enable_if<boost::is_scalar<typename GridType::ValueType>,
void>::type
4579 const GridType& grid,
4580 std::vector<Vec3s>& points,
4581 std::vector<Vec3I>& triangles,
4582 std::vector<Vec4I>& quads,
4595 tbb::parallel_for(tbb::blocked_range<size_t>(0, points.size()), ptnCpy);
4602 size_t numQuads = 0, numTriangles = 0;
4604 openvdb::tools::PolygonPool& polygons = polygonPoolList[n];
4605 numTriangles += polygons.numTriangles();
4606 numQuads += polygons.numQuads();
4610 triangles.resize(numTriangles);
4612 quads.resize(numQuads);
4616 size_t qIdx = 0, tIdx = 0;
4618 openvdb::tools::PolygonPool& polygons = polygonPoolList[n];
4620 for (
size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4621 quads[qIdx++] = polygons.quad(i);
4624 for (
size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
4625 triangles[tIdx++] = polygons.triangle(i);
4631 template<
typename Gr
idType>
4632 inline typename boost::disable_if<boost::is_scalar<typename GridType::ValueType>,
void>::type
4635 std::vector<Vec3s>&,
4636 std::vector<Vec3I>&,
4637 std::vector<Vec4I>&,
4645 template<
typename Gr
idType>
4648 const GridType& grid,
4649 std::vector<Vec3s>& points,
4650 std::vector<Vec3I>& triangles,
4651 std::vector<Vec4I>& quads,
4655 doVolumeToMesh(grid, points, triangles, quads, isovalue, adaptivity);
4659 template<
typename Gr
idType>
4662 const GridType& grid,
4663 std::vector<Vec3s>& points,
4664 std::vector<Vec4I>& quads,
4667 std::vector<Vec3I> triangles;
4679 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
void addLeaf(LeafNodeT *leaf)
Add the specified leaf to this tree, possibly creating a child branch in the process. If the leaf node already exists, replace it.
Definition: ValueAccessor.h:328
LeafNodeT * probeLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains voxel (x, y, z), or NULL if no such node exists...
Definition: ValueAccessor.h:378
OPENVDB_API const Index32 INVALID_IDX
boost::shared_ptr< const TreeBase > ConstPtr
Definition: Tree.h:67
Abstract base class for typed grids.
Definition: Grid.h:103
int16_t Int16
Definition: Types.h:60
Vec4< int32_t > Vec4i
Definition: Vec4.h:543
T & z()
Definition: Vec3.h:99
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
Gradient operators defined in index space of various orders.
Definition: Operators.h:122
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:109
math::Transform & transform()
Return a reference to this grid's transform, which might be shared with other grids.
Definition: Grid.h:335
int getValueDepth(const Coord &xyz) const
Definition: ValueAccessor.h:229
boost::shared_ptr< TreeBase > Ptr
Definition: Tree.h:66
bool diagonalizeSymmetricMatrix(const Mat3< T > &input, Mat3< T > &Q, Vec3< T > &D, unsigned int MAX_ITERATIONS=250)
Use Jacobi iterations to decompose a symmetric 3x3 matrix (diagonalize and compute eigenvectors) ...
Definition: Mat3.h:757
void setValueOn(const Coord &xyz, const ValueType &value)
Set the value of the voxel at the given coordinates and mark the voxel as active. ...
Definition: ValueAccessor.h:246
Defined various multi-threaded utility functions for trees.
size_t leafCount() const
Return the number of leaf nodes.
Definition: LeafManager.h:299
void foreach(const LeafOp &op, bool threaded=true, size_t grainSize=1)
Threaded method that applies a user-supplied functor to each leaf node in the LeafManager.
Definition: LeafManager.h:476
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: ValueAccessor.h:210
#define OPENVDB_VERSION_NAME
Definition: version.h:43
Index32 Index
Definition: Types.h:59
Definition: Exceptions.h:87
const LeafNodeT * probeConstLeaf(const Coord &xyz) const
Return a pointer to the leaf node that contains voxel (x, y, z), or NULL if no such node exists...
Definition: ValueAccessor.h:383
T dot(const Vec3< T > &v) const
Dot product.
Definition: Vec3.h:203
bool normalize(T eps=T(1.0e-7))
this = normalized this
Definition: Vec3.h:348
Definition: Exceptions.h:39
OPENVDB_API Hermite min(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
Vec3< double > Vec3d
Definition: Vec3.h:629
Mat3 transpose() const
returns transpose of this
Definition: Mat3.h:455
Mat3< double > Mat3d
Definition: Mat3.h:666
LeafNodeT * touchLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains voxel (x, y, z). If no such node exists, create one, but preserve the values and active states of all voxels.
Definition: ValueAccessor.h:347
Coord nearestCoord(const Vec3d &voxelCoord)
Return voxelCoord rounded to the closest integer coordinates.
Definition: Util.h:56
OPENVDB_API Hermite max(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:54
uint32_t Index32
Definition: Types.h:57
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
Vec3< float > Vec3s
Definition: Vec3.h:628
bool probeValue(const Coord &xyz, ValueType &value) const
Return the active state of the voxel as well as its value.
Definition: ValueAccessor.h:220
Base class for typed trees.
Definition: Tree.h:63
boost::shared_ptr< const GridBase > ConstPtr
Definition: Grid.h:107