38 #ifndef OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED 39 #define OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED 41 #include <openvdb/Platform.h> 42 #include <openvdb/math/Operators.h> 43 #include <openvdb/tree/ValueAccessor.h> 44 #include <openvdb/util/Util.h> 46 #include <boost/scoped_array.hpp> 47 #include <boost/scoped_ptr.hpp> 48 #include <boost/type_traits/is_scalar.hpp> 49 #include <boost/utility/enable_if.hpp> 50 #include <boost/math/special_functions/fpclassify.hpp> 52 #include <tbb/blocked_range.h> 53 #include <tbb/parallel_for.h> 54 #include <tbb/parallel_reduce.h> 55 #include <tbb/task_scheduler_init.h> 83 template<
typename Gr
idType>
87 std::vector<Vec3s>& points,
88 std::vector<Vec4I>& quads,
89 double isovalue = 0.0);
104 template<
typename Gr
idType>
107 const GridType& grid,
108 std::vector<Vec3s>& points,
109 std::vector<Vec3I>& triangles,
110 std::vector<Vec4I>& quads,
111 double isovalue = 0.0,
112 double adaptivity = 0.0,
129 inline PolygonPool(
const size_t numQuads,
const size_t numTriangles);
133 inline void resetQuads(
size_t size);
134 inline void clearQuads();
136 inline void resetTriangles(
size_t size);
137 inline void clearTriangles();
142 const size_t&
numQuads()
const {
return mNumQuads; }
157 const char&
quadFlags(
size_t n)
const {
return mQuadFlags[n]; }
166 inline bool trimQuads(
const size_t n,
bool reallocate =
false);
167 inline bool trimTrinagles(
const size_t n,
bool reallocate =
false);
173 size_t mNumQuads, mNumTriangles;
174 boost::scoped_array<openvdb::Vec4I> mQuads;
175 boost::scoped_array<openvdb::Vec3I> mTriangles;
176 boost::scoped_array<char> mQuadFlags, mTriangleFlags;
204 const size_t& pointListSize()
const;
205 PointList& pointList();
207 const size_t& polygonPoolListSize()
const;
208 PolygonPoolList& polygonPoolList();
209 const PolygonPoolList& polygonPoolList()
const;
211 std::vector<uint8_t>& pointFlags();
212 const std::vector<uint8_t>& pointFlags()
const;
220 template<
typename InputGr
idType>
221 void operator()(
const InputGridType&);
272 PolygonPoolList mPolygons;
274 size_t mPointListSize, mSeamPointListSize, mPolygonPoolListSize;
275 double mIsovalue, mPrimAdaptivity, mSecAdaptivity;
282 bool mInvertSurfaceMask, mRelaxDisorientedTriangles;
284 boost::scoped_array<uint32_t> mQuantizedSeamPoints;
285 std::vector<uint8_t> mPointFlags;
299 const std::vector<Vec3d>& points,
300 const std::vector<Vec3d>& normals)
306 if (points.empty())
return avgPos;
308 for (
size_t n = 0, N = points.size(); n < N; ++n) {
312 avgPos /= double(points.size());
316 double m00=0,m01=0,m02=0,
323 for (
size_t n = 0, N = points.size(); n < N; ++n) {
325 const Vec3d& n_ref = normals[n];
328 m00 += n_ref[0] * n_ref[0];
329 m11 += n_ref[1] * n_ref[1];
330 m22 += n_ref[2] * n_ref[2];
332 m01 += n_ref[0] * n_ref[1];
333 m02 += n_ref[0] * n_ref[2];
334 m12 += n_ref[1] * n_ref[2];
337 rhs += n_ref * n_ref.
dot(points[n] - avgPos);
362 Mat3d D = Mat3d::identity();
365 double tolerance =
std::max(std::abs(eigenValues[0]), std::abs(eigenValues[1]));
366 tolerance =
std::max(tolerance, std::abs(eigenValues[2]));
370 for (
int i = 0; i < 3; ++i ) {
371 if (std::abs(eigenValues[i]) < tolerance) {
375 D[i][i] = 1.0 / eigenValues[i];
381 Mat3d pseudoInv = eigenVectors * D * eigenVectors.
transpose();
382 return avgPos + pseudoInv * rhs;
396 namespace volume_to_mesh_internal {
402 #ifdef OPENVDB_HAS_CXX11 403 typedef std::unique_ptr<T>
type;
410 template<
typename ValueType>
413 FillArray(ValueType* array,
const ValueType& v) : mArray(array), mValue(v) { }
415 void operator()(
const tbb::blocked_range<size_t>& range)
const {
416 const ValueType v = mValue;
417 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
427 template<
typename ValueType>
429 fillArray(ValueType* array,
const ValueType& val,
const size_t length)
431 const size_t grainSize = length / tbb::task_scheduler_init::default_num_threads();
432 const tbb::blocked_range<size_t> range(0, length, grainSize);
444 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,
445 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,
446 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,
447 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,
448 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,
449 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,
450 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,
451 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};
456 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,
457 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,
458 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,
459 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,
460 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,
461 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,
462 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,
463 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};
470 {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},
471 {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},
472 {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},
473 {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},
474 {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},
475 {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},
476 {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},
477 {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},
478 {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},
479 {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},
480 {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},
481 {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},
482 {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},
483 {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},
484 {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},
485 {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},
486 {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},
487 {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},
488 {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},
489 {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},
490 {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},
491 {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},
492 {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},
493 {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},
494 {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},
495 {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},
496 {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},
497 {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},
498 {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},
499 {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},
500 {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},
501 {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},
502 {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},
503 {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},
504 {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},
505 {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},
506 {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},
507 {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},
508 {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},
509 {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},
510 {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},
511 {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},
512 {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},
513 {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},
514 {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},
515 {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},
516 {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},
517 {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},
518 {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},
519 {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},
520 {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},
521 {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},
522 {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},
523 {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},
524 {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},
525 {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},
526 {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},
527 {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},
528 {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},
529 {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},
530 {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},
531 {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},
532 {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},
533 {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},
534 {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},
535 {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},
536 {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},
537 {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},
538 {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},
539 {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},
540 {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},
541 {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},
542 {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},
543 {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},
544 {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},
545 {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},
546 {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},
547 {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},
548 {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},
549 {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},
550 {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},
551 {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},
552 {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},
553 {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},
554 {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},
555 {0,0,0,0,0,0,0,0,0,0,0,0,0}};
564 double epsilon = 0.001)
567 Vec3d normal = (p2-p0).cross(p1-p3);
569 const Vec3d centroid = (p0 + p1 + p2 + p3);
570 const double d = centroid.
dot(normal) * 0.25;
574 double absDist = std::abs(p0.dot(normal) - d);
575 if (absDist > epsilon)
return false;
577 absDist = std::abs(p1.dot(normal) - d);
578 if (absDist > epsilon)
return false;
580 absDist = std::abs(p2.dot(normal) - d);
581 if (absDist > epsilon)
return false;
583 absDist = std::abs(p3.dot(normal) - d);
584 if (absDist > epsilon)
return false;
608 assert(!(v.x() > 1.0) && !(v.y() > 1.0) && !(v.z() > 1.0));
609 assert(!(v.x() < 0.0) && !(v.y() < 0.0) && !(v.z() < 0.0));
624 v.y() = double(data & MASK_FIRST_10_BITS) * 0.0009775171;
626 v.x() = double(data & MASK_FIRST_10_BITS) * 0.0009775171;
650 template<
typename AccessorT>
655 values[0] = accessor.getValue(ijk);
657 values[1] = accessor.getValue(ijk);
659 values[2] = accessor.getValue(ijk);
661 values[3] = accessor.getValue(ijk);
663 values[4] = accessor.getValue(ijk);
665 values[5] = accessor.getValue(ijk);
667 values[6] = accessor.getValue(ijk);
669 values[7] = accessor.getValue(ijk);
673 template<
typename LeafT>
678 values[0] = leaf.getValue(offset);
679 values[3] = leaf.getValue(offset + 1);
680 values[4] = leaf.getValue(offset + LeafT::DIM);
681 values[7] = leaf.getValue(offset + LeafT::DIM + 1);
682 values[1] = leaf.getValue(offset + (LeafT::DIM * LeafT::DIM));
683 values[2] = leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + 1);
684 values[5] = leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM);
685 values[6] = leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM + 1);
689 template<
typename ValueType>
702 return uint8_t(signs);
708 template<
typename AccessorT>
710 evalCellSigns(
const AccessorT& accessor,
const Coord& ijk,
typename AccessorT::ValueType iso)
714 if (
isInsideValue(accessor.getValue(coord), iso)) signs |= 1u;
716 if (
isInsideValue(accessor.getValue(coord), iso)) signs |= 2u;
718 if (
isInsideValue(accessor.getValue(coord), iso)) signs |= 4u;
720 if (
isInsideValue(accessor.getValue(coord), iso)) signs |= 8u;
721 coord[1] += 1; coord[2] = ijk[2];
722 if (
isInsideValue(accessor.getValue(coord), iso)) signs |= 16u;
724 if (
isInsideValue(accessor.getValue(coord), iso)) signs |= 32u;
726 if (
isInsideValue(accessor.getValue(coord), iso)) signs |= 64u;
728 if (
isInsideValue(accessor.getValue(coord), iso)) signs |= 128u;
729 return uint8_t(signs);
735 template<
typename LeafT>
745 if (
isInsideValue(leaf.getValue(offset + 1), iso)) signs |= 8u;
748 if (
isInsideValue(leaf.getValue(offset + LeafT::DIM), iso)) signs |= 16u;
751 if (
isInsideValue(leaf.getValue(offset + LeafT::DIM + 1), iso)) signs |= 128u;
754 if (
isInsideValue(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) ), iso)) signs |= 2u;
757 if (
isInsideValue(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + 1), iso)) signs |= 4u;
760 if (
isInsideValue(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM), iso)) signs |= 32u;
763 if (
isInsideValue(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM + 1), iso)) signs |= 64u;
765 return uint8_t(signs);
771 template<
class AccessorT>
774 const AccessorT& acc, Coord ijk,
typename AccessorT::ValueType iso)
779 if (sAmbiguousFace[
evalCellSigns(acc, ijk, iso)] == 3) signs = uint8_t(~signs);
783 if (sAmbiguousFace[
evalCellSigns(acc, ijk, iso)] == 4) signs = uint8_t(~signs);
787 if (sAmbiguousFace[
evalCellSigns(acc, ijk, iso)] == 1) signs = uint8_t(~signs);
791 if (sAmbiguousFace[
evalCellSigns(acc, ijk, iso)] == 2) signs = uint8_t(~signs);
795 if (sAmbiguousFace[
evalCellSigns(acc, ijk, iso)] == 6) signs = uint8_t(~signs);
799 if (sAmbiguousFace[
evalCellSigns(acc, ijk, iso)] == 5) signs = uint8_t(~signs);
807 template<
class AccessorT>
810 typename AccessorT::ValueType isovalue,
const int dim)
823 coord[1] += dim; coord[2] = ijk[2];
834 if (p[0]) signs |= 1u;
835 if (p[1]) signs |= 2u;
836 if (p[2]) signs |= 4u;
837 if (p[3]) signs |= 8u;
838 if (p[4]) signs |= 16u;
839 if (p[5]) signs |= 32u;
840 if (p[6]) signs |= 64u;
841 if (p[7]) signs |= 128u;
842 if (!sAdaptable[signs])
return true;
847 int i = ijk[0], ip = ijk[0] + hDim, ipp = ijk[0] + dim;
848 int j = ijk[1], jp = ijk[1] + hDim, jpp = ijk[1] + dim;
849 int k = ijk[2], kp = ijk[2] + hDim, kpp = ijk[2] + dim;
852 coord.reset(ip, j, k);
854 if (p[0] != m && p[1] != m)
return true;
857 coord.reset(ipp, j, kp);
859 if (p[1] != m && p[2] != m)
return true;
862 coord.reset(ip, j, kpp);
864 if (p[2] != m && p[3] != m)
return true;
867 coord.reset(i, j, kp);
869 if (p[0] != m && p[3] != m)
return true;
872 coord.reset(ip, jpp, k);
874 if (p[4] != m && p[5] != m)
return true;
877 coord.reset(ipp, jpp, kp);
879 if (p[5] != m && p[6] != m)
return true;
882 coord.reset(ip, jpp, kpp);
884 if (p[6] != m && p[7] != m)
return true;
887 coord.reset(i, jpp, kp);
889 if (p[7] != m && p[4] != m)
return true;
892 coord.reset(i, jp, k);
894 if (p[0] != m && p[4] != m)
return true;
897 coord.reset(ipp, jp, k);
899 if (p[1] != m && p[5] != m)
return true;
902 coord.reset(ipp, jp, kpp);
904 if (p[2] != m && p[6] != m)
return true;
908 coord.reset(i, jp, kpp);
910 if (p[3] != m && p[7] != m)
return true;
916 coord.reset(ip, jp, k);
918 if (p[0] != m && p[1] != m && p[4] != m && p[5] != m)
return true;
921 coord.reset(ipp, jp, kp);
923 if (p[1] != m && p[2] != m && p[5] != m && p[6] != m)
return true;
926 coord.reset(ip, jp, kpp);
928 if (p[2] != m && p[3] != m && p[6] != m && p[7] != m)
return true;
931 coord.reset(i, jp, kp);
933 if (p[0] != m && p[3] != m && p[4] != m && p[7] != m)
return true;
936 coord.reset(ip, j, kp);
938 if (p[0] != m && p[1] != m && p[2] != m && p[3] != m)
return true;
941 coord.reset(ip, jpp, kp);
943 if (p[4] != m && p[5] != m && p[6] != m && p[7] != m)
return true;
946 coord.reset(ip, jp, kp);
948 if (p[0] != m && p[1] != m && p[2] != m && p[3] != m &&
949 p[4] != m && p[5] != m && p[6] != m && p[7] != m)
return true;
958 template <
class LeafType>
960 mergeVoxels(LeafType& leaf,
const Coord& start,
int dim,
int regionId)
962 Coord ijk, end = start;
967 for (ijk[0] = start[0]; ijk[0] < end[0]; ++ijk[0]) {
968 for (ijk[1] = start[1]; ijk[1] < end[1]; ++ijk[1]) {
969 for (ijk[2] = start[2]; ijk[2] < end[2]; ++ijk[2]) {
970 leaf.setValueOnly(ijk, regionId);
979 template <
class LeafType>
982 typename LeafType::ValueType::value_type adaptivity)
984 if (adaptivity < 1e-6)
return false;
986 typedef typename LeafType::ValueType VecT;
987 Coord ijk, end = start;
992 std::vector<VecT> norms;
993 for (ijk[0] = start[0]; ijk[0] < end[0]; ++ijk[0]) {
994 for (ijk[1] = start[1]; ijk[1] < end[1]; ++ijk[1]) {
995 for (ijk[2] = start[2]; ijk[2] < end[2]; ++ijk[2]) {
997 if(!leaf.isValueOn(ijk))
continue;
998 norms.push_back(leaf.getValue(ijk));
1003 size_t N = norms.size();
1004 for (
size_t ni = 0; ni < N; ++ni) {
1005 VecT n_i = norms[ni];
1006 for (
size_t nj = 0; nj < N; ++nj) {
1007 VecT n_j = norms[nj];
1008 if ((1.0 - n_i.dot(n_j)) > adaptivity)
return false;
1019 inline double evalZeroCrossing(
double v0,
double v1,
double iso) {
return (iso - v0) / (v1 - v0); }
1023 template<
typename LeafT>
1027 values[0] = double(leaf.getValue(offset));
1028 values[3] = double(leaf.getValue(offset + 1));
1029 values[4] = double(leaf.getValue(offset + LeafT::DIM));
1030 values[7] = double(leaf.getValue(offset + LeafT::DIM + 1));
1031 values[1] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM)));
1032 values[2] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + 1));
1033 values[5] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM));
1034 values[6] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM + 1));
1039 template<
typename AccessorT>
1044 values[0] = double(acc.getValue(coord));
1047 values[1] = double(acc.getValue(coord));
1050 values[2] = double(acc.getValue(coord));
1053 values[3] = double(acc.getValue(coord));
1055 coord[1] += 1; coord[2] = ijk[2];
1056 values[4] = double(acc.getValue(coord));
1059 values[5] = double(acc.getValue(coord));
1062 values[6] = double(acc.getValue(coord));
1065 values[7] = double(acc.getValue(coord));
1072 unsigned char edgeGroup,
double iso)
1074 Vec3d avg(0.0, 0.0, 0.0);
1077 if (sEdgeGroupTable[signs][1] == edgeGroup) {
1082 if (sEdgeGroupTable[signs][2] == edgeGroup) {
1088 if (sEdgeGroupTable[signs][3] == edgeGroup) {
1094 if (sEdgeGroupTable[signs][4] == edgeGroup) {
1099 if (sEdgeGroupTable[signs][5] == edgeGroup) {
1105 if (sEdgeGroupTable[signs][6] == edgeGroup) {
1112 if (sEdgeGroupTable[signs][7] == edgeGroup) {
1119 if (sEdgeGroupTable[signs][8] == edgeGroup) {
1125 if (sEdgeGroupTable[signs][9] == edgeGroup) {
1130 if (sEdgeGroupTable[signs][10] == edgeGroup) {
1136 if (sEdgeGroupTable[signs][11] == edgeGroup) {
1143 if (sEdgeGroupTable[signs][12] == edgeGroup) {
1150 double w = 1.0 / double(samples);
1164 unsigned char signsMask,
unsigned char edgeGroup,
double iso)
1166 avg =
Vec3d(0.0, 0.0, 0.0);
1169 if (sEdgeGroupTable[signs][1] == edgeGroup
1170 && sEdgeGroupTable[signsMask][1] == 0) {
1175 if (sEdgeGroupTable[signs][2] == edgeGroup
1176 && sEdgeGroupTable[signsMask][2] == 0) {
1182 if (sEdgeGroupTable[signs][3] == edgeGroup
1183 && sEdgeGroupTable[signsMask][3] == 0) {
1189 if (sEdgeGroupTable[signs][4] == edgeGroup
1190 && sEdgeGroupTable[signsMask][4] == 0) {
1195 if (sEdgeGroupTable[signs][5] == edgeGroup
1196 && sEdgeGroupTable[signsMask][5] == 0) {
1202 if (sEdgeGroupTable[signs][6] == edgeGroup
1203 && sEdgeGroupTable[signsMask][6] == 0) {
1210 if (sEdgeGroupTable[signs][7] == edgeGroup
1211 && sEdgeGroupTable[signsMask][7] == 0) {
1218 if (sEdgeGroupTable[signs][8] == edgeGroup
1219 && sEdgeGroupTable[signsMask][8] == 0) {
1225 if (sEdgeGroupTable[signs][9] == edgeGroup
1226 && sEdgeGroupTable[signsMask][9] == 0) {
1231 if (sEdgeGroupTable[signs][10] == edgeGroup
1232 && sEdgeGroupTable[signsMask][10] == 0) {
1238 if (sEdgeGroupTable[signs][11] == edgeGroup
1239 && sEdgeGroupTable[signsMask][11] == 0) {
1246 if (sEdgeGroupTable[signs][12] == edgeGroup
1247 && sEdgeGroupTable[signsMask][12] == 0) {
1254 double w = 1.0 / double(samples);
1268 unsigned char signs,
unsigned char edgeGroup,
double iso)
1270 std::vector<Vec3d> samples;
1273 std::vector<double> weights;
1276 Vec3d avg(0.0, 0.0, 0.0);
1278 if (sEdgeGroupTable[signs][1] == edgeGroup) {
1283 samples.push_back(avg);
1284 weights.push_back((avg-p).lengthSqr());
1287 if (sEdgeGroupTable[signs][2] == edgeGroup) {
1292 samples.push_back(avg);
1293 weights.push_back((avg-p).lengthSqr());
1296 if (sEdgeGroupTable[signs][3] == edgeGroup) {
1301 samples.push_back(avg);
1302 weights.push_back((avg-p).lengthSqr());
1305 if (sEdgeGroupTable[signs][4] == edgeGroup) {
1310 samples.push_back(avg);
1311 weights.push_back((avg-p).lengthSqr());
1314 if (sEdgeGroupTable[signs][5] == edgeGroup) {
1319 samples.push_back(avg);
1320 weights.push_back((avg-p).lengthSqr());
1323 if (sEdgeGroupTable[signs][6] == edgeGroup) {
1328 samples.push_back(avg);
1329 weights.push_back((avg-p).lengthSqr());
1332 if (sEdgeGroupTable[signs][7] == edgeGroup) {
1337 samples.push_back(avg);
1338 weights.push_back((avg-p).lengthSqr());
1341 if (sEdgeGroupTable[signs][8] == edgeGroup) {
1346 samples.push_back(avg);
1347 weights.push_back((avg-p).lengthSqr());
1350 if (sEdgeGroupTable[signs][9] == edgeGroup) {
1355 samples.push_back(avg);
1356 weights.push_back((avg-p).lengthSqr());
1359 if (sEdgeGroupTable[signs][10] == edgeGroup) {
1364 samples.push_back(avg);
1365 weights.push_back((avg-p).lengthSqr());
1368 if (sEdgeGroupTable[signs][11] == edgeGroup) {
1373 samples.push_back(avg);
1374 weights.push_back((avg-p).lengthSqr());
1377 if (sEdgeGroupTable[signs][12] == edgeGroup) {
1382 samples.push_back(avg);
1383 weights.push_back((avg-p).lengthSqr());
1390 for (
size_t i = 0, I = weights.size(); i < I; ++i) {
1391 minWeight =
std::min(minWeight, weights[i]);
1392 maxWeight =
std::max(maxWeight, weights[i]);
1395 const double offset = maxWeight + minWeight * 0.1;
1396 for (
size_t i = 0, I = weights.size(); i < I; ++i) {
1397 weights[i] = offset - weights[i];
1401 double weightSum = 0.0;
1402 for (
size_t i = 0, I = weights.size(); i < I; ++i) {
1403 weightSum += weights[i];
1410 if (samples.size() > 1) {
1411 for (
size_t i = 0, I = samples.size(); i < I; ++i) {
1412 avg += samples[i] * (weights[i] / weightSum);
1415 avg = samples.front();
1426 const std::vector<double>& values,
unsigned char signs,
double iso)
1428 for (
size_t n = 1, N = sEdgeGroupTable[signs][0] + 1; n < N; ++n) {
1429 points.push_back(
computePoint(values, signs, uint8_t(n), iso));
1438 matchEdgeGroup(
unsigned char groupId,
unsigned char lhsSigns,
unsigned char rhsSigns)
1441 for (
size_t i = 1; i <= 12; ++i) {
1442 if (sEdgeGroupTable[lhsSigns][i] == groupId && sEdgeGroupTable[rhsSigns][i] != 0) {
1443 id = sEdgeGroupTable[rhsSigns][i];
1457 const std::vector<double>& lhsValues,
const std::vector<double>& rhsValues,
1458 unsigned char lhsSigns,
unsigned char rhsSigns,
1459 double iso,
size_t pointIdx,
const uint32_t * seamPointArray)
1461 for (
size_t n = 1, N = sEdgeGroupTable[lhsSigns][0] + 1; n < N; ++n) {
1467 const unsigned char e = uint8_t(
id);
1468 const uint32_t& quantizedPoint = seamPointArray[pointIdx + (
id - 1)];
1473 weightedPointMask.push_back(
true);
1475 points.push_back(
computePoint(rhsValues, rhsSigns, e, iso));
1476 weightedPointMask.push_back(
false);
1480 points.push_back(
computePoint(lhsValues, lhsSigns, uint8_t(n), iso));
1481 weightedPointMask.push_back(
false);
1487 template <
typename InputTreeType>
1493 typedef typename InputTreeType::template ValueConverter<Int16>::Type
Int16TreeType;
1500 const InputTreeType& inputTree,
1501 const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
1502 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
1503 const boost::scoped_array<Index32>& leafNodeOffsets,
1507 void setRefData(
const InputTreeType& refInputTree,
1508 const Index32TreeType& refPointIndexTree,
1509 const Int16TreeType& refSignFlagsTree,
1510 const uint32_t * quantizedSeamLinePoints,
1511 uint8_t * seamLinePointsFlags);
1513 void operator()(
const tbb::blocked_range<size_t>&)
const;
1516 Vec3s *
const mPoints;
1517 InputTreeType
const *
const mInputTree;
1518 Index32LeafNodeType *
const *
const mPointIndexNodes;
1519 Int16LeafNodeType
const *
const *
const mSignFlagsNodes;
1520 Index32 const *
const mNodeOffsets;
1522 double const mIsovalue;
1524 InputTreeType
const * mRefInputTree;
1525 Index32TreeType
const * mRefPointIndexTree;
1526 Int16TreeType
const * mRefSignFlagsTree;
1527 uint32_t
const * mQuantizedSeamLinePoints;
1528 uint8_t * mSeamLinePointsFlags;
1532 template <
typename InputTreeType>
1535 const InputTreeType& inputTree,
1536 const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
1537 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
1538 const boost::scoped_array<Index32>& leafNodeOffsets,
1541 : mPoints(pointArray)
1542 , mInputTree(&inputTree)
1543 , mPointIndexNodes(pointIndexLeafNodes.empty() ? NULL : &pointIndexLeafNodes.front())
1544 , mSignFlagsNodes(signFlagsLeafNodes.empty() ? NULL : &signFlagsLeafNodes.front())
1545 , mNodeOffsets(leafNodeOffsets.get())
1548 , mRefInputTree(NULL)
1549 , mRefPointIndexTree(NULL)
1550 , mRefSignFlagsTree(NULL)
1551 , mQuantizedSeamLinePoints(NULL)
1552 , mSeamLinePointsFlags(NULL)
1556 template <
typename InputTreeType>
1559 const InputTreeType& refInputTree,
1562 const uint32_t * quantizedSeamLinePoints,
1563 uint8_t * seamLinePointsFlags)
1565 mRefInputTree = &refInputTree;
1566 mRefPointIndexTree = &refPointIndexTree;
1567 mRefSignFlagsTree = &refSignFlagsTree;
1568 mQuantizedSeamLinePoints = quantizedSeamLinePoints;
1569 mSeamLinePointsFlags = seamLinePointsFlags;
1572 template <
typename InputTreeType>
1580 typedef typename Index32TreeType::ValueType IndexType;
1582 typedef std::vector<Index> IndexArray;
1583 typedef std::map<IndexType, IndexArray> IndexArrayMap;
1585 InputTreeAccessor inputAcc(*mInputTree);
1589 std::vector<Vec3d> points(4);
1590 std::vector<bool> weightedPointMask(4);
1591 std::vector<double> values(8), refValues(8);
1592 const double iso = mIsovalue;
1596 boost::scoped_ptr<InputTreeAccessor> refInputAcc;
1597 boost::scoped_ptr<Index32TreeAccessor> refPointIndexAcc;
1598 boost::scoped_ptr<Int16TreeAccessor> refSignFlagsAcc;
1600 const bool hasReferenceData = mRefInputTree && mRefPointIndexTree && mRefSignFlagsTree;
1602 if (hasReferenceData) {
1603 refInputAcc.reset(
new InputTreeAccessor(*mRefInputTree));
1604 refPointIndexAcc.reset(
new Index32TreeAccessor(*mRefPointIndexTree));
1605 refSignFlagsAcc.reset(
new Int16TreeAccessor(*mRefSignFlagsTree));
1608 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
1611 const Coord& origin = pointIndexNode.origin();
1621 if (hasReferenceData) {
1622 refInputNode = refInputAcc->probeConstLeaf(origin);
1623 refPointIndexNode = refPointIndexAcc->probeConstLeaf(origin);
1624 refSignFlagsNode = refSignFlagsAcc->probeConstLeaf(origin);
1627 IndexType pointOffset = IndexType(mNodeOffsets[n]);
1628 IndexArrayMap regions;
1630 for (
typename Index32LeafNodeType::ValueOnIter it = pointIndexNode.beginValueOn(); it; ++it) {
1632 const Index offset = it.pos();
1634 const IndexType
id = it.getValue();
1637 regions[id].push_back(offset);
1642 pointIndexNode.setValueOnly(offset, pointOffset);
1644 const Int16 flags = signFlagsNode.getValue(offset);
1645 uint8_t signs = uint8_t(
SIGNS & flags);
1646 uint8_t refSigns = 0;
1648 if ((flags &
SEAM) && refPointIndexNode && refSignFlagsNode) {
1649 if (refSignFlagsNode->isValueOn(offset)) {
1650 refSigns = uint8_t(
SIGNS & refSignFlagsNode->getValue(offset));
1654 ijk = Index32LeafNodeType::offsetToLocalCoord(offset);
1656 const bool inclusiveCell = inputNode &&
1657 ijk[0] < int(Index32LeafNodeType::DIM - 1) &&
1658 ijk[1] < int(Index32LeafNodeType::DIM - 1) &&
1659 ijk[2] < int(Index32LeafNodeType::DIM - 1);
1667 weightedPointMask.clear();
1669 if (refSigns == 0) {
1677 computeCellPoints(points, weightedPointMask, values, refValues, signs, refSigns,
1678 iso, refPointIndexNode->getValue(offset), mQuantizedSeamLinePoints);
1681 xyz[0] = double(ijk[0]);
1682 xyz[1] = double(ijk[1]);
1683 xyz[2] = double(ijk[2]);
1685 for (
size_t i = 0, I = points.size(); i < I; ++i) {
1687 Vec3d& point = points[i];
1690 if (!boost::math::isfinite(point[0]) ||
1691 !boost::math::isfinite(point[1]) ||
1692 !boost::math::isfinite(point[2])) {
1694 " Hint: Check the input and consider using the \"Diagnostics\" tool " 1695 "to detect and resolve the NaNs.");
1701 Vec3s& pos = mPoints[pointOffset];
1702 pos[0] = float(point[0]);
1703 pos[1] = float(point[1]);
1704 pos[2] = float(point[2]);
1706 if (mSeamLinePointsFlags && !weightedPointMask.empty() && weightedPointMask[i]) {
1707 mSeamLinePointsFlags[pointOffset] = uint8_t(1);
1715 for (
typename IndexArrayMap::iterator it = regions.begin(); it != regions.end(); ++it) {
1717 Vec3d avg(0.0), point(0.0);
1720 const IndexArray& voxels = it->second;
1721 for (
size_t i = 0, I = voxels.size(); i < I; ++i) {
1723 const Index offset = voxels[i];
1724 ijk = Index32LeafNodeType::offsetToLocalCoord(offset);
1726 const bool inclusiveCell = inputNode &&
1727 ijk[0] < int(Index32LeafNodeType::DIM - 1) &&
1728 ijk[1] < int(Index32LeafNodeType::DIM - 1) &&
1729 ijk[2] < int(Index32LeafNodeType::DIM - 1);
1733 pointIndexNode.setValueOnly(offset, pointOffset);
1735 uint8_t signs = uint8_t(
SIGNS & signFlagsNode.getValue(offset));
1743 avg[0] += double(ijk[0]) + points[0][0];
1744 avg[1] += double(ijk[1]) + points[0][1];
1745 avg[2] += double(ijk[2]) + points[0][2];
1751 double w = 1.0 / double(count);
1759 Vec3s& pos = mPoints[pointOffset];
1760 pos[0] = float(avg[0]);
1761 pos[1] = float(avg[1]);
1762 pos[2] = float(avg[2]);
1773 template <
typename InputTreeType>
1779 typedef typename InputTreeType::template ValueConverter<Int16>::Type
Int16TreeType;
1786 const InputTreeType& inputTree,
1787 const Index32TreeType& refPointIndexTree,
1788 const Int16TreeType& refSignFlagsTree,
1789 uint32_t * quantizedPoints,
1791 : mSignFlagsNodes(signFlagsLeafNodes.empty() ? NULL : &signFlagsLeafNodes.front())
1792 , mInputTree(&inputTree)
1793 , mRefPointIndexTree(&refPointIndexTree)
1794 , mRefSignFlagsTree(&refSignFlagsTree)
1795 , mQuantizedPoints(quantizedPoints)
1806 std::vector<double> values(8);
1807 const double iso = double(mIsovalue);
1811 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
1813 const Int16LeafNodeType& signFlagsNode = *mSignFlagsNodes[n];
1814 const Coord& origin = signFlagsNode.origin();
1816 const Int16LeafNodeType * refSignNode = signFlagsTreeAcc.
probeConstLeaf(origin);
1817 if (!refSignNode)
continue;
1819 const Index32LeafNodeType * refPointIndexNode = pointIndexTreeAcc.
probeConstLeaf(origin);
1820 if (!refPointIndexNode)
continue;
1822 const InputLeafNodeType * inputNode = inputTreeAcc.
probeConstLeaf(origin);
1824 for (
typename Int16LeafNodeType::ValueOnCIter it = signFlagsNode.cbeginValueOn(); it; ++it) {
1826 const Index offset = it.pos();
1828 ijk = Index32LeafNodeType::offsetToLocalCoord(offset);
1830 const bool inclusiveCell = inputNode &&
1831 ijk[0] < int(Index32LeafNodeType::DIM - 1) &&
1832 ijk[1] < int(Index32LeafNodeType::DIM - 1) &&
1833 ijk[2] < int(Index32LeafNodeType::DIM - 1);
1837 if ((it.getValue() &
SEAM) && refSignNode->isValueOn(offset)) {
1839 uint8_t lhsSigns = uint8_t(
SIGNS & it.getValue());
1840 uint8_t rhsSigns = uint8_t(
SIGNS & refSignNode->getValue(offset));
1843 if (inclusiveCell) {
1850 for (
unsigned i = 1, I = sEdgeGroupTable[lhsSigns][0] + 1; i < I; ++i) {
1856 uint32_t& data = mQuantizedPoints[refPointIndexNode->getValue(offset) + (
id - 1)];
1861 pos, values, lhsSigns, rhsSigns, uint8_t(i), iso);
1878 Int16LeafNodeType
const *
const *
const mSignFlagsNodes;
1879 InputTreeType
const *
const mInputTree;
1880 Index32TreeType
const *
const mRefPointIndexTree;
1881 Int16TreeType
const *
const mRefSignFlagsTree;
1882 uint32_t *
const mQuantizedPoints;
1883 InputValueType
const mIsovalue;
1887 template <
typename TreeType>
1893 const TreeType& refSignFlagsTree)
1894 : mSignFlagsNodes(signFlagsLeafNodes.empty() ? NULL : &signFlagsLeafNodes.front())
1895 , mRefSignFlagsTree(&refSignFlagsTree)
1903 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
1905 LeafNodeType& signFlagsNode = *mSignFlagsNodes[n];
1906 const Coord& origin = signFlagsNode.origin();
1908 const LeafNodeType * refSignNode = refSignFlagsTreeAcc.
probeConstLeaf(origin);
1909 if (!refSignNode)
continue;
1911 for (
typename LeafNodeType::ValueOnCIter it = signFlagsNode.cbeginValueOn(); it; ++it) {
1912 const Index offset = it.pos();
1914 uint8_t rhsSigns = uint8_t(refSignNode->getValue(offset) &
SIGNS);
1916 if (sEdgeGroupTable[rhsSigns][0] > 0) {
1918 const typename LeafNodeType::ValueType value = it.getValue();
1919 uint8_t lhsSigns = uint8_t(value &
SIGNS);
1921 if (rhsSigns != lhsSigns) {
1922 signFlagsNode.setValueOnly(offset, value |
SEAM);
1932 LeafNodeType *
const *
const mSignFlagsNodes;
1933 TreeType
const *
const mRefSignFlagsTree;
1937 template <
typename BoolTreeType,
typename SignDataType>
1946 const BoolTreeType& maskTree)
1947 : mSignFlagsNodes(signFlagsLeafNodes.empty() ? NULL : &signFlagsLeafNodes.front())
1948 , mMaskTree(&maskTree)
1956 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
1958 SignDataLeafNodeType& signFlagsNode = *mSignFlagsNodes[n];
1959 const Coord& origin = signFlagsNode.origin();
1961 const BoolLeafNodeType * maskNode = maskAcc.
probeConstLeaf(origin);
1962 if (!maskNode)
continue;
1964 typedef typename SignDataLeafNodeType::ValueOnCIter ValueOnCIter;
1966 for (ValueOnCIter it = signFlagsNode.cbeginValueOn(); it; ++it) {
1967 const Index offset = it.pos();
1969 if (maskNode->isValueOn(offset)) {
1970 signFlagsNode.setValueOnly(offset, it.getValue() |
SEAM);
1979 SignDataLeafNodeType *
const *
const mSignFlagsNodes;
1980 BoolTreeType
const *
const mMaskTree;
1984 template <
typename TreeType>
1988 typedef typename TreeType::template ValueConverter<bool>::Type
BoolTreeType;
1991 const TreeType& signFlagsTree,
1993 : mSignFlagsNodes(signFlagsLeafNodes.empty() ? NULL : &signFlagsLeafNodes.front())
1994 , mSignFlagsTree(&signFlagsTree)
2001 : mSignFlagsNodes(rhs.mSignFlagsNodes)
2002 , mSignFlagsTree(rhs.mSignFlagsTree)
2012 typedef typename LeafNodeType::ValueOnCIter ValueOnCIter;
2013 typedef typename LeafNodeType::ValueType ValueType;
2019 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
2021 LeafNodeType& signFlagsNode = *mSignFlagsNodes[n];
2024 for (ValueOnCIter it = signFlagsNode.cbeginValueOn(); it; ++it) {
2026 const ValueType flags = it.getValue();
2028 if (!(flags &
SEAM) && (flags &
EDGES)) {
2030 ijk = it.getCoord();
2032 bool isSeamLineVoxel =
false;
2034 if (flags &
XEDGE) {
2038 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) &
SEAM);
2040 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) &
SEAM);
2044 if (!isSeamLineVoxel && flags &
YEDGE) {
2046 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) &
SEAM);
2048 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) &
SEAM);
2050 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) &
SEAM);
2054 if (!isSeamLineVoxel && flags &
ZEDGE) {
2056 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) &
SEAM);
2058 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) &
SEAM);
2060 isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.
getValue(ijk) &
SEAM);
2064 if (isSeamLineVoxel) {
2065 maskAcc.
setValue(it.getCoord(),
true);
2074 LeafNodeType *
const *
const mSignFlagsNodes;
2075 TreeType
const *
const mSignFlagsTree;
2076 BoolTreeType mTempMask;
2077 BoolTreeType *
const mMask;
2081 template<
typename SignDataTreeType>
2085 typedef typename SignDataTreeType::ValueType SignDataType;
2086 typedef typename SignDataTreeType::LeafNodeType SignDataLeafNodeType;
2087 typedef typename SignDataTreeType::template ValueConverter<bool>::Type BoolTreeType;
2089 std::vector<SignDataLeafNodeType*> signFlagsLeafNodes;
2090 signFlagsTree.getNodes(signFlagsLeafNodes);
2092 const tbb::blocked_range<size_t> nodeRange(0, signFlagsLeafNodes.size());
2094 tbb::parallel_for(nodeRange,
2097 BoolTreeType seamLineMaskTree(
false);
2100 maskSeamLine(signFlagsLeafNodes, signFlagsTree, seamLineMaskTree);
2102 tbb::parallel_reduce(nodeRange, maskSeamLine);
2104 tbb::parallel_for(nodeRange,
2112 template <
typename InputGr
idType>
2119 typedef typename InputTreeType::template ValueConverter<float>::Type
FloatTreeType;
2123 typedef typename InputTreeType::template ValueConverter<Int16>::Type
Int16TreeType;
2129 typedef typename InputTreeType::template ValueConverter<bool>::Type
BoolTreeType;
2133 const Index32TreeType& pointIndexTree,
2134 const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
2135 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
2138 bool invertSurfaceOrientation);
2142 mSpatialAdaptivityTree = &grid.
tree();
2143 mSpatialAdaptivityTransform = &grid.
transform();
2153 mRefSignFlagsTree = &signFlagsData;
2154 mInternalAdaptivity = internalAdaptivity;
2157 void operator()(
const tbb::blocked_range<size_t>&)
const;
2160 InputTreeType
const *
const mInputTree;
2163 Index32TreeType
const *
const mPointIndexTree;
2164 Index32LeafNodeType *
const *
const mPointIndexNodes;
2165 Int16LeafNodeType
const *
const *
const mSignFlagsNodes;
2167 InputValueType mIsovalue;
2168 float mSurfaceAdaptivity, mInternalAdaptivity;
2169 bool mInvertSurfaceOrientation;
2171 FloatTreeType
const * mSpatialAdaptivityTree;
2172 BoolTreeType
const * mMaskTree;
2173 Int16TreeType
const * mRefSignFlagsTree;
2178 template <
typename InputGr
idType>
2180 const InputGridType& inputGrid,
2182 const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
2183 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
2186 bool invertSurfaceOrientation)
2187 : mInputTree(&inputGrid.tree())
2188 , mInputTransform(&inputGrid.transform())
2189 , mPointIndexTree(&pointIndexTree)
2190 , mPointIndexNodes(pointIndexLeafNodes.empty() ? NULL : &pointIndexLeafNodes.front())
2191 , mSignFlagsNodes(signFlagsLeafNodes.empty() ? NULL : &signFlagsLeafNodes.front())
2193 , mSurfaceAdaptivity(adaptivity)
2194 , mInternalAdaptivity(adaptivity)
2195 , mInvertSurfaceOrientation(invertSurfaceOrientation)
2196 , mSpatialAdaptivityTree(NULL)
2198 , mRefSignFlagsTree(NULL)
2199 , mSpatialAdaptivityTransform(NULL)
2204 template <
typename InputGr
idType>
2209 typedef typename InputLeafNodeType::template ValueConverter<Vec3sType>::Type Vec3sLeafNodeType;
2217 boost::scoped_ptr<FloatTreeAccessor> spatialAdaptivityAcc;
2218 if (mSpatialAdaptivityTree && mSpatialAdaptivityTransform) {
2219 spatialAdaptivityAcc.reset(
new FloatTreeAccessor(*mSpatialAdaptivityTree));
2222 boost::scoped_ptr<BoolTreeAccessor> maskAcc;
2224 maskAcc.reset(
new BoolTreeAccessor(*mMaskTree));
2227 boost::scoped_ptr<Int16TreeAccessor> refSignFlagsAcc;
2228 if (mRefSignFlagsTree) {
2229 refSignFlagsAcc.reset(
new Int16TreeAccessor(*mRefSignFlagsTree));
2232 InputTreeAccessor inputAcc(*mInputTree);
2233 Index32TreeAccessor pointIndexAcc(*mPointIndexTree);
2237 const bool invertGradientDir = mInvertSurfaceOrientation || isBoolValue<InputValueType>();
2238 boost::scoped_ptr<Vec3sLeafNodeType> gradientNode;
2241 const int LeafDim = InputLeafNodeType::DIM;
2243 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
2245 mask.setValuesOff();
2250 const Coord& origin = pointIndexNode.origin();
2252 end[0] = origin[0] + LeafDim;
2253 end[1] = origin[1] + LeafDim;
2254 end[2] = origin[2] + LeafDim;
2259 if (maskLeaf != NULL) {
2260 for (
typename BoolLeafNodeType::ValueOnCIter it = maskLeaf->cbeginValueOn(); it; ++it) {
2261 mask.setActiveState(it.getCoord() & ~1u,
true);
2266 float adaptivity = (refSignFlagsAcc && !refSignFlagsAcc->probeConstLeaf(origin)) ?
2267 mInternalAdaptivity : mSurfaceAdaptivity;
2269 bool useGradients = adaptivity < 1.0f;
2274 if (spatialAdaptivityAcc) {
2275 useGradients =
false;
2276 for (
Index offset = 0; offset < FloatLeafNodeType::NUM_VALUES; ++offset) {
2277 ijk = adaptivityLeaf.offsetToGlobalCoord(offset);
2279 float weight = spatialAdaptivityAcc->getValue(ijk);
2280 float adaptivityValue = weight * adaptivity;
2281 if (adaptivityValue < 1.0f) useGradients =
true;
2282 adaptivityLeaf.setValueOnly(offset, adaptivityValue);
2287 for (
typename Int16LeafNodeType::ValueOnCIter it = signFlagsNode.cbeginValueOn(); it; ++it) {
2289 const Int16 flags = it.getValue();
2290 const unsigned char signs =
static_cast<unsigned char>(
SIGNS & int(flags));
2292 if ((flags &
SEAM) || !sAdaptable[signs] || sEdgeGroupTable[signs][0] > 1) {
2294 mask.setActiveState(it.getCoord() & ~1u,
true);
2296 }
else if (flags &
EDGES) {
2298 bool maskRegion =
false;
2300 ijk = it.getCoord();
2301 if (!pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2303 if (!maskRegion && flags &
XEDGE) {
2305 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2307 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2309 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2313 if (!maskRegion && flags &
YEDGE) {
2315 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2317 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2319 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2323 if (!maskRegion && flags &
ZEDGE) {
2325 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2327 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2329 if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion =
true;
2334 mask.setActiveState(it.getCoord() & ~1u,
true);
2341 for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
2342 for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
2343 for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
2344 if (!mask.isValueOn(ijk) &
isNonManifold(inputAcc, ijk, mIsovalue, dim)) {
2345 mask.setActiveState(ijk,
true);
2356 gradientNode->setValuesOff();
2358 gradientNode.
reset(
new Vec3sLeafNodeType());
2361 for (
typename Int16LeafNodeType::ValueOnCIter it = signFlagsNode.cbeginValueOn(); it; ++it) {
2362 ijk = it.getCoord();
2363 if (!mask.isValueOn(ijk & ~1u)) {
2367 if (invertGradientDir) {
2371 gradientNode->setValueOn(it.pos(), dir);
2378 for ( ; dim <= LeafDim; dim = dim << 1) {
2379 const unsigned coordMask = ~((dim << 1) - 1);
2380 for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
2381 for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
2382 for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
2384 adaptivity = adaptivityLeaf.getValue(ijk);
2386 if (mask.isValueOn(ijk) ||
isNonManifold(inputAcc, ijk, mIsovalue, dim)
2387 || (useGradients && !
isMergable(*gradientNode, ijk, dim, adaptivity)) ) {
2388 mask.setActiveState(ijk & coordMask,
true);
2390 mergeVoxels(pointIndexNode, ijk, dim, regionId++);
2410 mPolygonPool = &quadPool;
2415 template<
typename IndexType>
2419 mPolygonPool->quad(mIdx) = verts;
2421 Vec4I& quad = mPolygonPool->quad(mIdx);
2427 mPolygonPool->quadFlags(mIdx) = flags;
2433 mPolygonPool->trimQuads(mIdx);
2449 mPolygonPool = &polygonPool;
2451 mPolygonPool->resetTriangles(upperBound);
2457 template<
typename IndexType>
2460 if (verts[0] != verts[1] && verts[0] != verts[2] && verts[0] != verts[3]
2461 && verts[1] != verts[2] && verts[1] != verts[3] && verts[2] != verts[3]) {
2462 mPolygonPool->quadFlags(mQuadIdx) = flags;
2463 addQuad(verts, reverse);
2465 verts[0] == verts[3] &&
2466 verts[1] != verts[2] &&
2467 verts[1] != verts[0] &&
2468 verts[2] != verts[0]) {
2469 mPolygonPool->triangleFlags(mTriangleIdx) = flags;
2470 addTriangle(verts[0], verts[1], verts[2], reverse);
2472 verts[1] == verts[2] &&
2473 verts[0] != verts[3] &&
2474 verts[0] != verts[1] &&
2475 verts[3] != verts[1]) {
2476 mPolygonPool->triangleFlags(mTriangleIdx) = flags;
2477 addTriangle(verts[0], verts[1], verts[3], reverse);
2479 verts[0] == verts[1] &&
2480 verts[2] != verts[3] &&
2481 verts[2] != verts[0] &&
2482 verts[3] != verts[0]) {
2483 mPolygonPool->triangleFlags(mTriangleIdx) = flags;
2484 addTriangle(verts[0], verts[2], verts[3], reverse);
2486 verts[2] == verts[3] &&
2487 verts[0] != verts[1] &&
2488 verts[0] != verts[2] &&
2489 verts[1] != verts[2]) {
2490 mPolygonPool->triangleFlags(mTriangleIdx) = flags;
2491 addTriangle(verts[0], verts[1], verts[2], reverse);
2498 mPolygonPool->trimQuads(mQuadIdx,
true);
2499 mPolygonPool->trimTrinagles(mTriangleIdx,
true);
2504 template<
typename IndexType>
2508 mPolygonPool->quad(mQuadIdx) = verts;
2510 Vec4I& quad = mPolygonPool->quad(mQuadIdx);
2519 void addTriangle(
unsigned v0,
unsigned v1,
unsigned v2,
bool reverse)
2521 Vec3I& prim = mPolygonPool->triangle(mTriangleIdx);
2535 size_t mQuadIdx, mTriangleIdx;
2540 template<
typename SignAccT,
typename IdxAccT,
typename PrimBuilder>
2543 const SignAccT& signAcc,
const IdxAccT& idxAcc, PrimBuilder& mesher)
2545 typedef typename IdxAccT::ValueType IndexType;
2548 const bool isActive = idxAcc.probeValue(ijk, v0);
2555 bool isInside = flags &
INSIDE;
2557 isInside = invertSurfaceOrientation ? !isInside : isInside;
2562 if (flags &
XEDGE) {
2564 quad[0] = v0 + offsets[0];
2568 bool activeValues = idxAcc.probeValue(coord, quad[1]);
2569 uint8_t cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2570 quad[1] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][5] - 1 : 0;
2574 activeValues = activeValues && idxAcc.probeValue(coord, quad[2]);
2575 cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2576 quad[2] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][7] - 1 : 0;
2580 activeValues = activeValues && idxAcc.probeValue(coord, quad[3]);
2581 cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2582 quad[3] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][3] - 1 : 0;
2585 mesher.addPrim(quad, isInside, tag[
bool(refFlags & XEDGE)]);
2592 if (flags &
YEDGE) {
2594 quad[0] = v0 + offsets[1];
2598 bool activeValues = idxAcc.probeValue(coord, quad[1]);
2599 uint8_t cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2600 quad[1] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][12] - 1 : 0;
2604 activeValues = activeValues && idxAcc.probeValue(coord, quad[2]);
2605 cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2606 quad[2] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][11] - 1 : 0;
2610 activeValues = activeValues && idxAcc.probeValue(coord, quad[3]);
2611 cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2612 quad[3] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][10] - 1 : 0;
2615 mesher.addPrim(quad, isInside, tag[
bool(refFlags & YEDGE)]);
2622 if (flags &
ZEDGE) {
2624 quad[0] = v0 + offsets[2];
2628 bool activeValues = idxAcc.probeValue(coord, quad[1]);
2629 uint8_t cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2630 quad[1] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][8] - 1 : 0;
2634 activeValues = activeValues && idxAcc.probeValue(coord, quad[2]);
2635 cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2636 quad[2] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][6] - 1 : 0;
2640 activeValues = activeValues && idxAcc.probeValue(coord, quad[3]);
2641 cell = uint8_t(
SIGNS & signAcc.getValue(coord));
2642 quad[3] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][2] - 1 : 0;
2645 mesher.addPrim(quad, !isInside, tag[
bool(refFlags & ZEDGE)]);
2654 template<
typename InputTreeType>
2658 typedef typename InputTreeType::template ValueConverter<bool>::Type
BoolTreeType;
2662 : mInputTree(&inputTree)
2666 , mTileArray(tileArray)
2671 : mInputTree(rhs.mInputTree)
2672 , mIsovalue(rhs.mIsovalue)
2675 , mTileArray(rhs.mTileArray)
2681 void operator()(
const tbb::blocked_range<size_t>&);
2685 InputValueType
const mIsovalue;
2686 BoolTreeType mTempMask;
2687 BoolTreeType *
const mMask;
2688 Vec4i const *
const mTileArray;
2692 template<
typename InputTreeType>
2698 CoordBBox region, bbox;
2701 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
2703 const Vec4i& tile = mTileArray[n];
2705 bbox.min()[0] = tile[0];
2706 bbox.min()[1] = tile[1];
2707 bbox.min()[2] = tile[2];
2708 bbox.max() = bbox.min();
2709 bbox.max().offset(tile[3]);
2714 const int valueDepth = inputTreeAcc.
getValueDepth(bbox.min());
2722 bool processRegion =
true;
2727 if (processRegion) {
2730 region.min()[0] = region.max()[0] = ijk[0];
2731 mMask->fill(region,
false);
2738 processRegion =
true;
2743 if (processRegion) {
2746 region.min()[0] = region.max()[0] = ijk[0];
2747 mMask->fill(region,
false);
2757 processRegion =
true;
2762 if (processRegion) {
2765 region.min()[1] = region.max()[1] = ijk[1];
2766 mMask->fill(region,
false);
2773 processRegion =
true;
2778 if (processRegion) {
2781 region.min()[1] = region.max()[1] = ijk[1];
2782 mMask->fill(region,
false);
2792 processRegion =
true;
2797 if (processRegion) {
2800 region.min()[2] = region.max()[2] = ijk[2];
2801 mMask->fill(region,
false);
2807 processRegion =
true;
2812 if (processRegion) {
2815 region.min()[2] = region.max()[2] = ijk[2];
2816 mMask->fill(region,
false);
2822 template<
typename InputTreeType>
2825 typename InputTreeType::template ValueConverter<bool>::Type& mask)
2827 typename InputTreeType::ValueOnCIter tileIter(inputTree);
2828 tileIter.setMaxDepth(InputTreeType::ValueOnCIter::LEAF_DEPTH - 1);
2830 size_t tileCount = 0;
2831 for ( ; tileIter; ++tileIter) {
2835 if (tileCount > 0) {
2836 boost::scoped_array<Vec4i> tiles(
new Vec4i[tileCount]);
2841 tileIter = inputTree.cbeginValueOn();
2842 tileIter.setMaxDepth(InputTreeType::ValueOnCIter::LEAF_DEPTH - 1);
2844 for (; tileIter; ++tileIter) {
2845 Vec4i& tile = tiles[index++];
2846 tileIter.getBoundingBox(bbox);
2847 tile[0] = bbox.min()[0];
2848 tile[1] = bbox.min()[1];
2849 tile[2] = bbox.min()[2];
2850 tile[3] = bbox.max()[0] - bbox.min()[0];
2854 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, tileCount), op);
2867 : mPointsIn(pointsIn) , mPointsOut(pointsOut)
2873 for (
size_t n = range.begin(); n < range.end(); ++n) {
2874 mPointsOut[n] = mPointsIn[n];
2879 const PointList& mPointsIn;
2880 std::vector<Vec3s>& mPointsOut;
2892 template<
typename LeafNodeType>
2893 void constructOffsetList();
2896 const IndexVector&
core()
const {
return mCore; }
2900 const IndexVector&
minX()
const {
return mMinX; }
2903 const IndexVector&
maxX()
const {
return mMaxX; }
2907 const IndexVector&
minY()
const {
return mMinY; }
2910 const IndexVector&
maxY()
const {
return mMaxY; }
2914 const IndexVector&
minZ()
const {
return mMinZ; }
2917 const IndexVector&
maxZ()
const {
return mMaxZ; }
2931 IndexVector mCore, mMinX, mMaxX, mMinY, mMaxY, mMinZ, mMaxZ,
2932 mInternalNeighborsX, mInternalNeighborsY, mInternalNeighborsZ;
2936 template<
typename LeafNodeType>
2942 mCore.reserve((LeafNodeType::DIM - 2) * (LeafNodeType::DIM - 2));
2944 for (
Index x = 1; x < (LeafNodeType::DIM - 1); ++x) {
2945 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2946 for (
Index y = 1; y < (LeafNodeType::DIM - 1); ++y) {
2947 const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
2948 for (
Index z = 1; z < (LeafNodeType::DIM - 1); ++z) {
2949 mCore.push_back(offsetXY + z);
2955 mInternalNeighborsX.clear();
2956 mInternalNeighborsX.reserve(LeafNodeType::SIZE - (LeafNodeType::DIM * LeafNodeType::DIM));
2958 for (
Index x = 0; x < (LeafNodeType::DIM - 1); ++x) {
2959 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2960 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
2961 const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
2962 for (
Index z = 0; z < LeafNodeType::DIM; ++z) {
2963 mInternalNeighborsX.push_back(offsetXY + z);
2969 mInternalNeighborsY.clear();
2970 mInternalNeighborsY.reserve(LeafNodeType::SIZE - (LeafNodeType::DIM * LeafNodeType::DIM));
2972 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
2973 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2974 for (
Index y = 0; y < (LeafNodeType::DIM - 1); ++y) {
2975 const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
2976 for (
Index z = 0; z < LeafNodeType::DIM; ++z) {
2977 mInternalNeighborsY.push_back(offsetXY + z);
2983 mInternalNeighborsZ.clear();
2984 mInternalNeighborsZ.reserve(LeafNodeType::SIZE - (LeafNodeType::DIM * LeafNodeType::DIM));
2986 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
2987 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2988 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
2989 const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
2990 for (
Index z = 0; z < (LeafNodeType::DIM - 1); ++z) {
2991 mInternalNeighborsZ.push_back(offsetXY + z);
2998 mMinX.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
3000 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
3001 const Index offsetXY = (y << LeafNodeType::LOG2DIM);
3002 for (
Index z = 0; z < LeafNodeType::DIM; ++z) {
3003 mMinX.push_back(offsetXY + z);
3010 mMaxX.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
3012 const Index offsetX = (LeafNodeType::DIM - 1) << (2 * LeafNodeType::LOG2DIM);
3013 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
3014 const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
3015 for (
Index z = 0; z < LeafNodeType::DIM; ++z) {
3016 mMaxX.push_back(offsetXY + z);
3023 mMinY.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
3025 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
3026 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
3027 for (
Index z = 0; z < (LeafNodeType::DIM - 1); ++z) {
3028 mMinY.push_back(offsetX + z);
3035 mMaxY.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
3037 const Index offsetY = (LeafNodeType::DIM - 1) << LeafNodeType::LOG2DIM;
3038 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
3039 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
3040 for (
Index z = 0; z < (LeafNodeType::DIM - 1); ++z) {
3041 mMaxY.push_back(offsetX + offsetY + z);
3048 mMinZ.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
3050 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
3051 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
3052 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
3053 const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
3054 mMinZ.push_back(offsetXY);
3061 mMaxZ.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
3063 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
3064 const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
3065 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
3066 const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
3067 mMaxZ.push_back(offsetXY + (LeafNodeType::DIM - 1));
3078 template<
typename AccessorT,
int _AXIS>
3081 enum { AXIS = _AXIS };
3086 void set(Coord ijk) {
3088 acc.setActiveState(ijk);
3090 acc.setActiveState(ijk);
3092 acc.setActiveState(ijk);
3094 acc.setActiveState(ijk);
3095 }
else if (_AXIS == 1) {
3096 acc.setActiveState(ijk);
3098 acc.setActiveState(ijk);
3100 acc.setActiveState(ijk);
3102 acc.setActiveState(ijk);
3104 acc.setActiveState(ijk);
3106 acc.setActiveState(ijk);
3108 acc.setActiveState(ijk);
3110 acc.setActiveState(ijk);
3119 template<
typename VoxelEdgeAcc,
typename LeafNode>
3127 if (VoxelEdgeAcc::AXIS == 0) {
3128 nvo = LeafNode::DIM * LeafNode::DIM;
3130 }
else if (VoxelEdgeAcc::AXIS == 1) {
3131 nvo = LeafNode::DIM;
3135 for (
size_t n = 0, N = offsets->size(); n < N; ++n) {
3136 const Index& pos = (*offsets)[n];
3137 bool isActive = leafnode.isValueOn(pos) || leafnode.isValueOn(pos + nvo);
3138 if (isActive && (
isInsideValue(leafnode.getValue(pos), iso) !=
3140 edgeAcc.set(leafnode.offsetToGlobalCoord(pos));
3149 template<
typename LeafNode,
typename TreeAcc,
typename VoxelEdgeAcc>
3154 const std::vector<Index>* lhsOffsets = &voxels.
maxX();
3155 const std::vector<Index>* rhsOffsets = &voxels.
minX();
3156 Coord ijk = lhsNode.origin();
3158 if (VoxelEdgeAcc::AXIS == 0) {
3159 ijk[0] += LeafNode::DIM;
3160 }
else if (VoxelEdgeAcc::AXIS == 1) {
3161 ijk[1] += LeafNode::DIM;
3162 lhsOffsets = &voxels.
maxY();
3163 rhsOffsets = &voxels.
minY();
3164 }
else if (VoxelEdgeAcc::AXIS == 2) {
3165 ijk[2] += LeafNode::DIM;
3166 lhsOffsets = &voxels.
maxZ();
3167 rhsOffsets = &voxels.
minZ();
3170 typename LeafNode::ValueType value;
3171 const LeafNode* rhsNodePt = acc.probeConstLeaf(ijk);
3174 for (
size_t n = 0, N = lhsOffsets->size(); n < N; ++n) {
3175 const Index& pos = (*lhsOffsets)[n];
3176 bool isActive = lhsNode.isValueOn(pos) || rhsNodePt->isValueOn((*rhsOffsets)[n]);
3177 if (isActive && (
isInsideValue(lhsNode.getValue(pos), iso) !=
3178 isInsideValue(rhsNodePt->getValue((*rhsOffsets)[n]), iso))) {
3179 edgeAcc.set(lhsNode.offsetToGlobalCoord(pos));
3182 }
else if (!acc.probeValue(ijk, value)) {
3184 for (
size_t n = 0, N = lhsOffsets->size(); n < N; ++n) {
3185 const Index& pos = (*lhsOffsets)[n];
3186 if (lhsNode.isValueOn(pos) && (inside !=
isInsideValue(lhsNode.getValue(pos), iso))) {
3187 edgeAcc.set(lhsNode.offsetToGlobalCoord(pos));
3197 template<
typename LeafNode,
typename TreeAcc,
typename VoxelEdgeAcc>
3202 Coord ijk = leafnode.origin();
3203 if (VoxelEdgeAcc::AXIS == 0) --ijk[0];
3204 else if (VoxelEdgeAcc::AXIS == 1) --ijk[1];
3205 else if (VoxelEdgeAcc::AXIS == 2) --ijk[2];
3207 typename LeafNode::ValueType value;
3208 if (!acc.probeConstLeaf(ijk) && !acc.probeValue(ijk, value)) {
3215 for (
size_t n = 0, N = offsets->size(); n < N; ++n) {
3217 const Index& pos = (*offsets)[n];
3218 if (leafnode.isValueOn(pos) && (inside !=
isInsideValue(leafnode.getValue(pos), iso))) {
3220 ijk = leafnode.offsetToGlobalCoord(pos);
3221 if (VoxelEdgeAcc::AXIS == 0) --ijk[0];
3222 else if (VoxelEdgeAcc::AXIS == 1) --ijk[1];
3223 else if (VoxelEdgeAcc::AXIS == 2) --ijk[2];
3233 template<
typename InputTreeType>
3239 typedef typename InputTreeType::template ValueConverter<bool>::Type
BoolTreeType;
3243 const std::vector<const InputLeafNodeType*>& inputLeafNodes,
3244 BoolTreeType& intersectionTree,
3245 InputValueType iso);
3248 void operator()(
const tbb::blocked_range<size_t>&);
3250 mIntersectionAccessor.tree().merge(rhs.mIntersectionAccessor.
tree());
3255 InputLeafNodeType
const *
const *
const mInputNodes;
3257 BoolTreeType mIntersectionTree;
3263 InputValueType mIsovalue;
3267 template<
typename InputTreeType>
3270 const std::vector<const InputLeafNodeType*>& inputLeafNodes,
3273 : mInputAccessor(inputTree)
3274 , mInputNodes(inputLeafNodes.empty() ? NULL : &inputLeafNodes.front())
3275 , mIntersectionTree(false)
3276 , mIntersectionAccessor(intersectionTree)
3278 , mOffsets(&mOffsetData)
3285 template<
typename InputTreeType>
3288 : mInputAccessor(rhs.mInputAccessor.tree())
3289 , mInputNodes(rhs.mInputNodes)
3290 , mIntersectionTree(false)
3291 , mIntersectionAccessor(mIntersectionTree)
3293 , mOffsets(rhs.mOffsets)
3294 , mIsovalue(rhs.mIsovalue)
3299 template<
typename InputTreeType>
3307 for (
size_t n = range.begin(); n != range.end(); ++n) {
3338 template<
typename InputTreeType>
3341 typename InputTreeType::template ValueConverter<bool>::Type& intersectionTree,
3342 const InputTreeType& inputTree,
3343 typename InputTreeType::ValueType isovalue)
3346 typedef typename InputTreeType::template ValueConverter<bool>::Type
BoolTreeType;
3348 std::vector<const InputLeafNodeType*> inputLeafNodes;
3349 inputTree.getNodes(inputLeafNodes);
3352 inputTree, inputLeafNodes, intersectionTree, isovalue);
3354 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, inputLeafNodes.size()), op);
3363 template<
typename InputTreeType>
3369 typedef typename InputTreeType::template ValueConverter<bool>::Type
BoolTreeType;
3373 const InputTreeType& inputTree,
3374 const std::vector<BoolLeafNodeType*>& nodes,
3375 BoolTreeType& intersectionTree,
3376 InputValueType iso);
3379 void operator()(
const tbb::blocked_range<size_t>&);
3381 mIntersectionAccessor.tree().merge(rhs.mIntersectionAccessor.
tree());
3386 BoolLeafNodeType
const *
const *
const mNodes;
3388 BoolTreeType mIntersectionTree;
3391 InputValueType mIsovalue;
3395 template<
typename InputTreeType>
3397 const InputTreeType& inputTree,
3398 const std::vector<BoolLeafNodeType*>& nodes,
3401 : mInputAccessor(inputTree)
3402 , mNodes(nodes.empty() ? NULL : &nodes.front())
3403 , mIntersectionTree(false)
3404 , mIntersectionAccessor(intersectionTree)
3410 template<
typename InputTreeType>
3413 : mInputAccessor(rhs.mInputAccessor.tree())
3414 , mNodes(rhs.mNodes)
3415 , mIntersectionTree(false)
3416 , mIntersectionAccessor(mIntersectionTree)
3417 , mIsovalue(rhs.mIsovalue)
3422 template<
typename InputTreeType>
3433 for (
size_t n = range.begin(); n != range.end(); ++n) {
3437 for (
typename BoolLeafNodeType::ValueOnCIter it = node.cbeginValueOn(); it; ++it) {
3439 if (!it.getValue()) {
3441 ijk = it.getCoord();
3462 template<
typename BoolTreeType>
3468 const std::vector<BoolLeafNodeType*>& maskNodes,
3470 : mMaskTree(&maskTree)
3471 , mMaskNodes(maskNodes.empty() ? NULL : &maskNodes.front())
3472 , mTmpBorderTree(false)
3473 , mBorderTree(&borderTree)
3478 : mMaskTree(rhs.mMaskTree)
3479 , mMaskNodes(rhs.mMaskNodes)
3480 , mTmpBorderTree(false)
3481 , mBorderTree(&mTmpBorderTree)
3493 for (
size_t n = range.begin(); n != range.end(); ++n) {
3495 const BoolLeafNodeType& node = *mMaskNodes[n];
3497 for (
typename BoolLeafNodeType::ValueOnCIter it = node.cbeginValueOn(); it; ++it) {
3499 ijk = it.getCoord();
3501 const bool lhs = it.getValue();
3504 bool isEdgeVoxel =
false;
3507 isEdgeVoxel = (maskAcc.
probeValue(ijk, rhs) && lhs != rhs);
3510 isEdgeVoxel = isEdgeVoxel || (maskAcc.
probeValue(ijk, rhs) && lhs != rhs);
3513 isEdgeVoxel = isEdgeVoxel || (maskAcc.
probeValue(ijk, rhs) && lhs != rhs);
3516 isEdgeVoxel = isEdgeVoxel || (maskAcc.
probeValue(ijk, rhs) && lhs != rhs);
3520 isEdgeVoxel = isEdgeVoxel || (maskAcc.
probeValue(ijk, rhs) && lhs != rhs);
3523 isEdgeVoxel = isEdgeVoxel || (maskAcc.
probeValue(ijk, rhs) && lhs != rhs);
3526 isEdgeVoxel = isEdgeVoxel || (maskAcc.
probeValue(ijk, rhs) && lhs != rhs);
3538 BoolLeafNodeType
const *
const *
const mMaskNodes;
3545 template<
typename BoolTreeType>
3551 : mNodes(nodes.empty() ? NULL : &nodes.front())
3558 typedef typename BoolLeafNodeType::ValueOnIter ValueOnIter;
3562 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
3564 BoolLeafNodeType& node = *mNodes[n];
3566 const BoolLeafNodeType * maskNode = maskTreeAcc.
probeConstLeaf(node.origin());
3569 for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3570 const Index pos = it.pos();
3571 if (maskNode->getValue(pos)) {
3572 node.setValueOnly(pos,
true);
3580 BoolLeafNodeType *
const *
const mNodes;
3588 template<
typename BoolTreeType>
3595 : mNodes(nodes.empty() ? NULL : &nodes.front())
3597 , mInputTransform(inputTransform)
3598 , mMaskTransform(maskTransform)
3599 , mInvertMask(invert)
3605 typedef typename BoolLeafNodeType::ValueOnIter ValueOnIter;
3609 const bool matchingTransforms = mInputTransform == mMaskTransform;
3610 const bool maskState = mInvertMask;
3612 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
3614 BoolLeafNodeType& node = *mNodes[n];
3616 if (matchingTransforms) {
3618 const BoolLeafNodeType * maskNode = maskTreeAcc.
probeConstLeaf(node.origin());
3622 for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3623 const Index pos = it.pos();
3624 if (maskNode->isValueOn(pos) == maskState) {
3625 node.setValueOnly(pos,
true);
3631 if (maskTreeAcc.
isValueOn(node.origin()) == maskState) {
3632 for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3633 node.setValueOnly(it.pos(),
true);
3643 for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3645 ijk = mMaskTransform.worldToIndexCellCentered(
3646 mInputTransform.indexToWorld(it.getCoord()));
3648 if (maskTreeAcc.
isValueOn(ijk) == maskState) {
3649 node.setValueOnly(it.pos(),
true);
3658 BoolLeafNodeType *
const *
const mNodes;
3662 bool const mInvertMask;
3666 template<
typename InputGr
idType>
3669 typename InputGridType::TreeType::template ValueConverter<bool>::Type& intersectionTree,
3670 typename InputGridType::TreeType::template ValueConverter<bool>::Type& borderTree,
3671 const InputGridType& inputGrid,
3674 typename InputGridType::ValueType isovalue)
3676 typedef typename InputGridType::TreeType InputTreeType;
3679 typedef typename InputTreeType::template ValueConverter<bool>::Type
BoolTreeType;
3683 if (maskGrid && maskGrid->type() == BoolGridType::gridType()) {
3686 const InputTreeType& inputTree = inputGrid.tree();
3688 const BoolGridType * surfaceMask =
static_cast<const BoolGridType*
>(maskGrid.get());
3690 const BoolTreeType& maskTree = surfaceMask->tree();
3696 std::vector<BoolLeafNodeType*> intersectionLeafNodes;
3697 intersectionTree.getNodes(intersectionLeafNodes);
3699 tbb::parallel_for(tbb::blocked_range<size_t>(0, intersectionLeafNodes.size()),
3701 intersectionLeafNodes, maskTree, transform, maskTransform, invertMask));
3707 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, intersectionLeafNodes.size()), borderOp);
3712 BoolTreeType tmpIntersectionTree(
false);
3715 inputTree, intersectionLeafNodes, tmpIntersectionTree, isovalue);
3717 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, intersectionLeafNodes.size()), op);
3719 std::vector<BoolLeafNodeType*> tmpIntersectionLeafNodes;
3720 tmpIntersectionTree.getNodes(tmpIntersectionLeafNodes);
3722 tbb::parallel_for(tbb::blocked_range<size_t>(0, tmpIntersectionLeafNodes.size()),
3725 intersectionTree.clear();
3726 intersectionTree.merge(tmpIntersectionTree);
3734 template<
typename InputTreeType>
3742 typedef typename InputTreeType::template ValueConverter<Int16>::Type
Int16TreeType;
3747 const std::vector<const BoolLeafNodeType*>& intersectionLeafNodes,
3748 Int16TreeType& signFlagsTree,
3749 Index32TreeType& pointIndexTree,
3750 InputValueType iso);
3753 void operator()(
const tbb::blocked_range<size_t>&);
3755 mSignFlagsAccessor.tree().merge(rhs.mSignFlagsAccessor.
tree());
3756 mPointIndexAccessor.tree().merge(rhs.mPointIndexAccessor.
tree());
3761 BoolLeafNodeType
const *
const *
const mIntersectionNodes;
3763 Int16TreeType mSignFlagsTree;
3765 Index32TreeType mPointIndexTree;
3768 const InputValueType mIsovalue;
3772 template<
typename InputTreeType>
3774 const InputTreeType& inputTree,
3775 const std::vector<const BoolLeafNodeType*>& intersectionLeafNodes,
3779 : mInputAccessor(inputTree)
3780 , mIntersectionNodes(&intersectionLeafNodes.front())
3782 , mSignFlagsAccessor(signFlagsTree)
3783 , mPointIndexTree(boost::integer_traits<
Index32>::const_max)
3784 , mPointIndexAccessor(pointIndexTree)
3787 pointIndexTree.root().setBackground(boost::integer_traits<Index32>::const_max,
false);
3791 template<
typename InputTreeType>
3793 : mInputAccessor(rhs.mInputAccessor.tree())
3794 , mIntersectionNodes(rhs.mIntersectionNodes)
3796 , mSignFlagsAccessor(mSignFlagsTree)
3797 , mPointIndexTree(boost::integer_traits<
Index32>::const_max)
3798 , mPointIndexAccessor(mPointIndexTree)
3799 , mIsovalue(rhs.mIsovalue)
3804 template<
typename InputTreeType>
3808 typedef typename Int16TreeType::LeafNodeType Int16LeafNodeType;
3814 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
3817 const Coord& origin = maskNode.
origin();
3821 if (!signsNodePt.get()) signsNodePt.reset(
new Int16LeafNodeType(origin, 0));
3822 else signsNodePt->setOrigin(origin);
3824 bool updatedNode =
false;
3828 const Index pos = it.pos();
3842 if (signFlags != 0 && signFlags != 0xFF) {
3844 const bool inside = signFlags & 0x1;
3846 int edgeFlags = inside ?
INSIDE : 0;
3848 if (!it.getValue()) {
3849 edgeFlags |= inside != ((signFlags & 0x02) != 0) ?
XEDGE : 0;
3850 edgeFlags |= inside != ((signFlags & 0x10) != 0) ?
YEDGE : 0;
3851 edgeFlags |= inside != ((signFlags & 0x08) != 0) ?
ZEDGE : 0;
3854 const uint8_t ambiguousCellFlags = sAmbiguousFace[signFlags];
3855 if (ambiguousCellFlags != 0) {
3856 correctCellSigns(signFlags, ambiguousCellFlags, mInputAccessor, origin + ijk, mIsovalue);
3859 edgeFlags |= int(signFlags);
3861 signsNodePt->setValueOn(pos,
Int16(edgeFlags));
3867 typename Index32TreeType::LeafNodeType* idxNode = mPointIndexAccessor.
touchLeaf(origin);
3868 idxNode->topologyUnion(*signsNodePt);
3871 for (
typename Index32TreeType::LeafNodeType::ValueOnIter it = idxNode->beginValueOn(); it; ++it) {
3872 idxNode->setValueOnly(it.pos(), 0);
3875 mSignFlagsAccessor.
addLeaf(signsNodePt.release());
3881 template<
typename InputTreeType>
3884 typename InputTreeType::template ValueConverter<Int16>::Type& signFlagsTree,
3885 typename InputTreeType::template ValueConverter<Index32>::Type& pointIndexTree,
3886 const typename InputTreeType::template ValueConverter<bool>::Type& intersectionTree,
3887 const InputTreeType& inputTree,
3888 typename InputTreeType::ValueType isovalue)
3890 typedef typename InputTreeType::template ValueConverter<bool>::Type BoolTreeType;
3893 std::vector<const BoolLeafNodeType*> intersectionLeafNodes;
3894 intersectionTree.getNodes(intersectionLeafNodes);
3897 inputTree, intersectionLeafNodes, signFlagsTree, pointIndexTree, isovalue);
3899 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, intersectionLeafNodes.size()), op);
3906 template<Index32 LeafNodeLog2Dim>
3912 boost::scoped_array<Index32>& leafNodeCount)
3913 : mLeafNodes(leafNodes.empty() ? NULL : &leafNodes.front())
3914 , mData(leafNodeCount.get())
3920 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
3924 Int16 const * p = mLeafNodes[n]->buffer().data();
3925 Int16 const *
const endP = p + Int16LeafNodeType::SIZE;
3937 Int16LeafNodeType *
const *
const mLeafNodes;
3942 template<
typename Po
intIndexLeafNode>
3948 const std::vector<Int16LeafNodeType*>& signDataNodes,
3949 boost::scoped_array<Index32>& leafNodeCount)
3950 : mPointIndexNodes(pointIndexNodes.empty() ? NULL : &pointIndexNodes.front())
3951 , mSignDataNodes(signDataNodes.empty() ? NULL : &signDataNodes.front())
3952 , mData(leafNodeCount.get())
3960 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
3963 const Int16LeafNodeType& signNode = *mSignDataNodes[n];
3967 std::set<IndexType> uniqueRegions;
3974 count += size_t(sEdgeGroupTable[(
SIGNS & signNode.
getValue(it.pos()))][0]);
3976 uniqueRegions.insert(
id);
3980 mData[n] =
Index32(count + uniqueRegions.size());
3986 Int16LeafNodeType
const *
const *
const mSignDataNodes;
3991 template<
typename Po
intIndexLeafNode>
3996 MapPoints(
const std::vector<PointIndexLeafNode*>& pointIndexNodes,
3997 const std::vector<Int16LeafNodeType*>& signDataNodes,
3998 boost::scoped_array<Index32>& leafNodeCount)
3999 : mPointIndexNodes(pointIndexNodes.empty() ? NULL : &pointIndexNodes.front())
4000 , mSignDataNodes(signDataNodes.empty() ? NULL : &signDataNodes.front())
4001 , mData(leafNodeCount.get())
4007 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
4009 const Int16LeafNodeType& signNode = *mSignDataNodes[n];
4012 Index32 pointOffset = mData[n];
4015 const Index pos = it.pos();
4018 pointOffset +=
Index32(sEdgeGroupTable[signs][0]);
4025 Int16LeafNodeType
const *
const *
const mSignDataNodes;
4032 template<
typename TreeType,
typename PrimBuilder>
4043 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
4044 const Int16TreeType& signFlagsTree,
4045 const Index32TreeType& idxTree,
4046 PolygonPoolList& polygons,
4047 bool invertSurfaceOrientation);
4051 void operator()(
const tbb::blocked_range<size_t>&)
const;
4054 Int16LeafNodeType *
const *
const mSignFlagsLeafNodes;
4055 Int16TreeType
const *
const mSignFlagsTree;
4056 Int16TreeType
const * mRefSignFlagsTree;
4057 Index32TreeType
const *
const mIndexTree;
4058 PolygonPoolList *
const mPolygonPoolList;
4059 bool const mInvertSurfaceOrientation;
4063 template<
typename TreeType,
typename PrimBuilder>
4065 const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
4068 PolygonPoolList& polygons,
4069 bool invertSurfaceOrientation)
4070 : mSignFlagsLeafNodes(signFlagsLeafNodes.empty() ? NULL : &signFlagsLeafNodes.front())
4071 , mSignFlagsTree(&signFlagsTree)
4072 , mRefSignFlagsTree(NULL)
4073 , mIndexTree(&idxTree)
4074 , mPolygonPoolList(&polygons)
4075 , mInvertSurfaceOrientation(invertSurfaceOrientation)
4079 template<
typename InputTreeType,
typename PrimBuilder>
4084 Int16ValueAccessor signAcc(*mSignFlagsTree);
4088 const bool invertSurfaceOrientation = mInvertSurfaceOrientation;
4095 boost::scoped_ptr<Int16ValueAccessor> refSignAcc;
4096 if (mRefSignFlagsTree) refSignAcc.reset(
new Int16ValueAccessor(*mRefSignFlagsTree));
4098 for (
size_t n = range.begin(); n != range.end(); ++n) {
4101 origin = node.origin();
4105 typename Int16LeafNodeType::ValueOnCIter iter = node.cbeginValueOn();
4106 for (; iter; ++iter) {
4107 if (iter.getValue() &
XEDGE) ++edgeCount;
4108 if (iter.getValue() &
YEDGE) ++edgeCount;
4109 if (iter.getValue() &
ZEDGE) ++edgeCount;
4112 if(edgeCount == 0)
continue;
4114 mesher.init(edgeCount, (*mPolygonPoolList)[n]);
4119 if (!signleafPt || !idxLeafPt)
continue;
4123 if (refSignAcc) refSignLeafPt = refSignAcc->probeConstLeaf(origin);
4127 for (iter = node.cbeginValueOn(); iter; ++iter) {
4128 ijk = iter.getCoord();
4130 Int16 flags = iter.getValue();
4132 if (!(flags & 0xE00))
continue;
4135 if (refSignLeafPt) {
4136 refFlags = refSignLeafPt->getValue(iter.pos());
4143 const uint8_t cell = uint8_t(
SIGNS & flags);
4145 if (sEdgeGroupTable[cell][0] > 1) {
4146 offsets[0] = (sEdgeGroupTable[cell][1] - 1);
4147 offsets[1] = (sEdgeGroupTable[cell][9] - 1);
4148 offsets[2] = (sEdgeGroupTable[cell][4] - 1);
4151 if (ijk[0] > origin[0] && ijk[1] > origin[1] && ijk[2] > origin[2]) {
4153 flags, refFlags, offsets, ijk, *signleafPt, *idxLeafPt, mesher);
4156 flags, refFlags, offsets, ijk, signAcc, idxAcc, mesher);
4169 template<
typename T>
4172 CopyArray(T * outputArray,
const T * inputArray,
size_t outputOffset = 0)
4173 : mOutputArray(outputArray), mInputArray(inputArray), mOutputOffset(outputOffset)
4177 void operator()(
const tbb::blocked_range<size_t>& inputArrayRange)
const 4179 const size_t offset = mOutputOffset;
4180 for (
size_t n = inputArrayRange.begin(), N = inputArrayRange.end(); n < N; ++n) {
4181 mOutputArray[offset + n] = mInputArray[n];
4186 T *
const mOutputArray;
4187 T
const *
const mInputArray;
4188 size_t const mOutputOffset;
4195 const std::vector<uint8_t>& pointFlags,
4196 boost::scoped_array<openvdb::Vec3s>& points,
4197 boost::scoped_array<unsigned>& numQuadsToDivide)
4198 : mPolygonPoolList(&polygons)
4199 , mPointFlags(pointFlags.empty() ? NULL : &pointFlags.front())
4200 , mPoints(points.get())
4201 , mNumQuadsToDivide(numQuadsToDivide.get())
4207 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
4214 for (
size_t i = 0, I = polygons.
numQuads(); i < I; ++i) {
4222 const bool edgePoly = mPointFlags[quad[0]] || mPointFlags[quad[1]]
4223 || mPointFlags[quad[2]] || mPointFlags[quad[3]];
4225 if (!edgePoly)
continue;
4227 const Vec3s& p0 = mPoints[quad[0]];
4228 const Vec3s& p1 = mPoints[quad[1]];
4229 const Vec3s& p2 = mPoints[quad[2]];
4230 const Vec3s& p3 = mPoints[quad[3]];
4239 mNumQuadsToDivide[n] = count;
4244 PolygonPoolList *
const mPolygonPoolList;
4245 uint8_t
const *
const mPointFlags;
4246 Vec3s const *
const mPoints;
4247 unsigned *
const mNumQuadsToDivide;
4254 const boost::scoped_array<openvdb::Vec3s>& points,
4256 boost::scoped_array<openvdb::Vec3s>& centroids,
4257 boost::scoped_array<unsigned>& numQuadsToDivide,
4258 boost::scoped_array<unsigned>& centroidOffsets)
4259 : mPolygonPoolList(&polygons)
4260 , mPoints(points.get())
4261 , mCentroids(centroids.get())
4262 , mNumQuadsToDivide(numQuadsToDivide.get())
4263 , mCentroidOffsets(centroidOffsets.get())
4264 , mPointCount(pointCount)
4270 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
4274 const size_t nonplanarCount = size_t(mNumQuadsToDivide[n]);
4276 if (nonplanarCount > 0) {
4282 size_t offset = mCentroidOffsets[n];
4284 size_t triangleIdx = 0;
4286 for (
size_t i = 0, I = polygons.
numQuads(); i < I; ++i) {
4288 const char quadFlags = polygons.
quadFlags(i);
4291 unsigned newPointIdx = unsigned(offset + mPointCount);
4295 mCentroids[offset] = (mPoints[quad[0]] + mPoints[quad[1]] +
4296 mPoints[quad[2]] + mPoints[quad[3]]) * 0.25f;
4303 triangle[0] = quad[0];
4304 triangle[1] = newPointIdx;
4305 triangle[2] = quad[3];
4315 triangle[0] = quad[0];
4316 triangle[1] = quad[1];
4317 triangle[2] = newPointIdx;
4327 triangle[0] = quad[1];
4328 triangle[1] = quad[2];
4329 triangle[2] = newPointIdx;
4340 triangle[0] = quad[2];
4341 triangle[1] = quad[3];
4342 triangle[2] = newPointIdx;
4353 for (
size_t i = 0, I = polygons.
numTriangles(); i < I; ++i) {
4360 for (
size_t i = 0, I = polygons.
numQuads(); i < I; ++i) {
4364 tmpPolygons.
quad(quadIdx) = quad;
4370 polygons.
copy(tmpPolygons);
4376 PolygonPoolList *
const mPolygonPoolList;
4377 Vec3s const *
const mPoints;
4378 Vec3s *
const mCentroids;
4379 unsigned *
const mNumQuadsToDivide;
4380 unsigned *
const mCentroidOffsets;
4381 size_t const mPointCount;
4387 PolygonPoolList& polygonPoolList,
4388 size_t polygonPoolListSize,
4389 PointList& pointList,
4390 size_t& pointListSize,
4391 std::vector<uint8_t>& pointFlags)
4393 const tbb::blocked_range<size_t> polygonPoolListRange(0, polygonPoolListSize);
4395 boost::scoped_array<unsigned> numQuadsToDivide(
new unsigned[polygonPoolListSize]);
4397 tbb::parallel_for(polygonPoolListRange,
4400 boost::scoped_array<unsigned> centroidOffsets(
new unsigned[polygonPoolListSize]);
4402 size_t centroidCount = 0;
4406 for (
size_t n = 0, N = polygonPoolListSize; n < N; ++n) {
4407 centroidOffsets[n] = sum;
4408 sum += numQuadsToDivide[n];
4410 centroidCount = size_t(sum);
4413 boost::scoped_array<Vec3s> centroidList(
new Vec3s[centroidCount]);
4415 tbb::parallel_for(polygonPoolListRange,
4417 centroidList, numQuadsToDivide, centroidOffsets));
4419 if (centroidCount > 0) {
4421 const size_t newPointListSize = centroidCount + pointListSize;
4423 boost::scoped_array<openvdb::Vec3s> newPointList(
new openvdb::Vec3s[newPointListSize]);
4425 tbb::parallel_for(tbb::blocked_range<size_t>(0, pointListSize),
4428 tbb::parallel_for(tbb::blocked_range<size_t>(0, newPointListSize - pointListSize),
4431 pointListSize = newPointListSize;
4432 pointList.swap(newPointList);
4433 pointFlags.resize(pointListSize, 0);
4441 const std::vector<uint8_t>& pointFlags)
4442 : mPolygonPoolList(&polygons)
4443 , mPointFlags(pointFlags.empty() ? NULL : &pointFlags.front())
4449 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
4453 for (
size_t i = 0, I = polygons.
numQuads(); i < I; ++i) {
4461 const bool hasSeamLinePoint =
4462 mPointFlags[verts[0]] || mPointFlags[verts[1]] ||
4463 mPointFlags[verts[2]] || mPointFlags[verts[3]];
4465 if (!hasSeamLinePoint) {
4466 flags &= ~POLYFLAG_FRACTURE_SEAM;
4471 for (
size_t i = 0, I = polygons.
numTriangles(); i < I; ++i) {
4479 const bool hasSeamLinePoint =
4480 mPointFlags[verts[0]] || mPointFlags[verts[1]] || mPointFlags[verts[2]];
4482 if (!hasSeamLinePoint) {
4483 flags &= ~POLYFLAG_FRACTURE_SEAM;
4493 PolygonPoolList *
const mPolygonPoolList;
4494 uint8_t
const *
const mPointFlags;
4500 std::vector<uint8_t>& pointFlags)
4502 tbb::parallel_for(tbb::blocked_range<size_t>(0, polygonPoolListSize),
4510 template<
typename InputTreeType>
4514 const PointList& pointList, boost::scoped_array<uint8_t>& pointMask,
4516 : mInputTree(&inputTree)
4517 , mPolygonPoolList(&polygons)
4518 , mPointList(&pointList)
4519 , mPointMask(pointMask.get())
4520 , mTransform(transform)
4521 , mInvertSurfaceOrientation(invertSurfaceOrientation)
4527 typedef typename InputTreeType::LeafNodeType::ValueType ValueType;
4530 Vec3s centroid, normal;
4533 const bool invertGradientDir = mInvertSurfaceOrientation || isBoolValue<ValueType>();
4535 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
4537 const PolygonPool& polygons = (*mPolygonPoolList)[n];
4539 for (
size_t i = 0, I = polygons.
numTriangles(); i < I; ++i) {
4543 const Vec3s& v0 = (*mPointList)[verts[0]];
4544 const Vec3s& v1 = (*mPointList)[verts[1]];
4545 const Vec3s& v2 = (*mPointList)[verts[2]];
4547 normal = (v2 - v0).cross((v1 - v0));
4550 centroid = (v0 + v1 + v2) * (1.0f / 3.0f);
4551 ijk = mTransform.worldToIndexCellCentered(centroid);
4556 if (invertGradientDir) {
4561 if (dir.dot(normal) < -0.5f) {
4566 mPointMask[verts[0]] = 1;
4567 mPointMask[verts[1]] = 1;
4568 mPointMask[verts[2]] = 1;
4577 InputTreeType
const *
const mInputTree;
4578 PolygonPoolList
const *
const mPolygonPoolList;
4579 PointList
const *
const mPointList;
4580 uint8_t *
const mPointMask;
4582 bool const mInvertSurfaceOrientation;
4586 template<
typename InputTree>
4589 bool invertSurfaceOrientation,
4590 const InputTree& inputTree,
4592 PolygonPoolList& polygonPoolList,
4593 size_t polygonPoolListSize,
4594 PointList& pointList,
4595 const size_t pointListSize)
4597 const tbb::blocked_range<size_t> polygonPoolListRange(0, polygonPoolListSize);
4599 boost::scoped_array<uint8_t> pointMask(
new uint8_t[pointListSize]);
4600 fillArray(pointMask.get(), uint8_t(0), pointListSize);
4602 tbb::parallel_for(polygonPoolListRange,
4604 inputTree, polygonPoolList, pointList, pointMask, transform, invertSurfaceOrientation));
4606 boost::scoped_array<uint8_t> pointUpdates(
new uint8_t[pointListSize]);
4607 fillArray(pointUpdates.get(), uint8_t(0), pointListSize);
4609 boost::scoped_array<Vec3s> newPoints(
new Vec3s[pointListSize]);
4610 fillArray(newPoints.get(),
Vec3s(0.0f, 0.0f, 0.0f), pointListSize);
4612 for (
size_t n = 0, N = polygonPoolListSize; n < N; ++n) {
4616 for (
size_t i = 0, I = polygons.
numQuads(); i < I; ++i) {
4619 for (
int v = 0; v < 4; ++v) {
4621 const unsigned pointIdx = verts[v];
4623 if (pointMask[pointIdx] == 1) {
4625 newPoints[pointIdx] +=
4626 pointList[verts[0]] + pointList[verts[1]] +
4627 pointList[verts[2]] + pointList[verts[3]];
4629 pointUpdates[pointIdx] = uint8_t(pointUpdates[pointIdx] + 4);
4634 for (
size_t i = 0, I = polygons.
numTriangles(); i < I; ++i) {
4637 for (
int v = 0; v < 3; ++v) {
4639 const unsigned pointIdx = verts[v];
4641 if (pointMask[pointIdx] == 1) {
4642 newPoints[pointIdx] +=
4643 pointList[verts[0]] + pointList[verts[1]] + pointList[verts[2]];
4645 pointUpdates[pointIdx] = uint8_t(pointUpdates[pointIdx] + 3);
4651 for (
size_t n = 0, N = pointListSize; n < N; ++n) {
4652 if (pointUpdates[n] > 0) {
4653 const double weight = 1.0 / double(pointUpdates[n]);
4654 pointList[n] = newPoints[n] * float(weight);
4673 , mTriangleFlags(NULL)
4680 : mNumQuads(numQuads)
4681 , mNumTriangles(numTriangles)
4684 , mQuadFlags(new char[mNumQuads])
4685 , mTriangleFlags(new char[mNumTriangles])
4696 for (
size_t i = 0; i < mNumQuads; ++i) {
4697 mQuads[i] = rhs.mQuads[i];
4698 mQuadFlags[i] = rhs.mQuadFlags[i];
4701 for (
size_t i = 0; i < mNumTriangles; ++i) {
4702 mTriangles[i] = rhs.mTriangles[i];
4703 mTriangleFlags[i] = rhs.mTriangleFlags[i];
4713 mQuadFlags.reset(
new char[mNumQuads]);
4722 mQuadFlags.reset(NULL);
4729 mNumTriangles = size;
4731 mTriangleFlags.reset(
new char[mNumTriangles]);
4739 mTriangles.reset(NULL);
4740 mTriangleFlags.reset(NULL);
4747 if (!(n < mNumQuads))
return false;
4755 boost::scoped_array<openvdb::Vec4I> quads(
new openvdb::Vec4I[n]);
4756 boost::scoped_array<char> flags(
new char[n]);
4758 for (
size_t i = 0; i < n; ++i) {
4759 quads[i] = mQuads[i];
4760 flags[i] = mQuadFlags[i];
4764 mQuadFlags.swap(flags);
4776 if (!(n < mNumTriangles))
return false;
4781 mTriangles.reset(NULL);
4784 boost::scoped_array<openvdb::Vec3I> triangles(
new openvdb::Vec3I[n]);
4785 boost::scoped_array<char> flags(
new char[n]);
4787 for (
size_t i = 0; i < n; ++i) {
4788 triangles[i] = mTriangles[i];
4789 flags[i] = mTriangleFlags[i];
4792 mTriangles.swap(triangles);
4793 mTriangleFlags.swap(flags);
4809 , mSeamPointListSize(0)
4810 , mPolygonPoolListSize(0)
4811 , mIsovalue(isovalue)
4812 , mPrimAdaptivity(adaptivity)
4813 , mSecAdaptivity(0.0)
4815 , mSurfaceMaskGrid(
GridBase::ConstPtr())
4816 , mAdaptivityGrid(
GridBase::ConstPtr())
4817 , mAdaptivityMaskTree(
TreeBase::ConstPtr())
4820 , mInvertSurfaceMask(false)
4821 , mRelaxDisorientedTriangles(relaxDisorientedTriangles)
4822 , mQuantizedSeamPoints(NULL)
4835 inline const size_t&
4838 return mPointListSize;
4842 inline PolygonPoolList&
4849 inline const PolygonPoolList&
4856 inline const size_t&
4859 return mPolygonPoolListSize;
4867 mSecAdaptivity = secAdaptivity;
4872 mSeamPointListSize = 0;
4873 mQuantizedSeamPoints.reset(NULL);
4880 mSurfaceMaskGrid = mask;
4881 mInvertSurfaceMask = invertMask;
4888 mAdaptivityGrid = grid;
4895 mAdaptivityMaskTree = tree;
4899 inline std::vector<uint8_t>&
4906 inline const std::vector<uint8_t>&
4913 template<
typename InputGr
idType>
4919 typedef typename InputGridType::TreeType InputTreeType;
4920 typedef typename InputTreeType::LeafNodeType InputLeafNodeType;
4921 typedef typename InputLeafNodeType::ValueType InputValueType;
4926 typedef typename InputTreeType::template ValueConverter<float>::Type FloatTreeType;
4929 typedef typename InputTreeType::template ValueConverter<bool>::Type BoolTreeType;
4930 typedef typename BoolTreeType::LeafNodeType BoolLeafNodeType;
4933 typedef typename InputTreeType::template ValueConverter<Int16>::Type Int16TreeType;
4934 typedef typename Int16TreeType::LeafNodeType Int16LeafNodeType;
4936 typedef typename InputTreeType::template ValueConverter<Index32>::Type Index32TreeType;
4937 typedef typename Index32TreeType::LeafNodeType Index32LeafNodeType;
4944 mPolygonPoolListSize = 0;
4946 mPointFlags.clear();
4952 const InputValueType isovalue = InputValueType(mIsovalue);
4953 const float adaptivityThreshold = float(mPrimAdaptivity);
4954 const bool adaptive = mPrimAdaptivity > 1e-7 || mSecAdaptivity > 1e-7;
4960 const bool invertSurfaceOrientation = !volume_to_mesh_internal::isBoolValue<InputValueType>() &&
4966 const InputTreeType& inputTree = inputGrid.tree();
4968 BoolTreeType intersectionTree(
false), adaptivityMask(
false);
4970 if (mAdaptivityMaskTree && mAdaptivityMaskTree->type() == BoolTreeType::treeType()) {
4971 const BoolTreeType *refAdaptivityMask=
4972 static_cast<const BoolTreeType*
>(mAdaptivityMaskTree.get());
4973 adaptivityMask.topologyUnion(*refAdaptivityMask);
4976 Int16TreeType signFlagsTree(0);
4977 Index32TreeType pointIndexTree(boost::integer_traits<Index32>::const_max);
4983 intersectionTree, inputTree, isovalue);
4986 intersectionTree, adaptivityMask, inputGrid, mSurfaceMaskGrid, mInvertSurfaceMask, isovalue);
4988 if (intersectionTree.empty())
return;
4991 signFlagsTree, pointIndexTree, intersectionTree, inputTree, isovalue);
4993 intersectionTree.clear();
4995 std::vector<Index32LeafNodeType*> pointIndexLeafNodes;
4996 pointIndexTree.getNodes(pointIndexLeafNodes);
4998 std::vector<Int16LeafNodeType*> signFlagsLeafNodes;
4999 signFlagsTree.getNodes(signFlagsLeafNodes);
5001 const tbb::blocked_range<size_t> auxiliaryLeafNodeRange(0, signFlagsLeafNodes.size());
5006 Int16TreeType * refSignFlagsTree = NULL;
5007 Index32TreeType * refPointIndexTree = NULL;
5008 InputTreeType
const * refInputTree = NULL;
5010 if (mRefGrid && mRefGrid->type() == InputGridType::gridType()) {
5012 const InputGridType* refGrid =
static_cast<const InputGridType*
>(mRefGrid.get());
5013 refInputTree = &refGrid->tree();
5015 if (!mRefSignTree && !mRefIdxTree) {
5019 typename Int16TreeType::Ptr refSignFlagsTreePt(
new Int16TreeType(0));
5020 typename Index32TreeType::Ptr refPointIndexTreePt(
5021 new Index32TreeType(boost::integer_traits<Index32>::const_max));
5023 BoolTreeType refIntersectionTree(
false);
5026 refIntersectionTree, *refInputTree, isovalue);
5029 *refPointIndexTreePt, refIntersectionTree, *refInputTree, isovalue);
5031 mRefSignTree = refSignFlagsTreePt;
5032 mRefIdxTree = refPointIndexTreePt;
5035 if (mRefSignTree && mRefIdxTree) {
5039 refSignFlagsTree =
static_cast<Int16TreeType*
>(mRefSignTree.get());
5040 refPointIndexTree =
static_cast<Index32TreeType*
>(mRefIdxTree.get());
5044 if (refSignFlagsTree && refPointIndexTree) {
5050 if (mSeamPointListSize == 0) {
5054 std::vector<Int16LeafNodeType*> refSignFlagsLeafNodes;
5055 refSignFlagsTree->getNodes(refSignFlagsLeafNodes);
5057 boost::scoped_array<Index32> leafNodeOffsets(
new Index32[refSignFlagsLeafNodes.size()]);
5059 tbb::parallel_for(tbb::blocked_range<size_t>(0, refSignFlagsLeafNodes.size()),
5061 refSignFlagsLeafNodes, leafNodeOffsets));
5065 for (
size_t n = 0, N = refSignFlagsLeafNodes.size(); n < N; ++n) {
5066 const Index32 tmp = leafNodeOffsets[n];
5067 leafNodeOffsets[n] = count;
5070 mSeamPointListSize = size_t(count);
5073 if (mSeamPointListSize != 0) {
5075 mQuantizedSeamPoints.reset(
new uint32_t[mSeamPointListSize]);
5077 memset(mQuantizedSeamPoints.get(), 0,
sizeof(uint32_t) * mSeamPointListSize);
5079 std::vector<Index32LeafNodeType*> refPointIndexLeafNodes;
5080 refPointIndexTree->getNodes(refPointIndexLeafNodes);
5082 tbb::parallel_for(tbb::blocked_range<size_t>(0, refPointIndexLeafNodes.size()),
5084 refPointIndexLeafNodes, refSignFlagsLeafNodes, leafNodeOffsets));
5088 if (mSeamPointListSize != 0) {
5090 tbb::parallel_for(auxiliaryLeafNodeRange,
5092 signFlagsLeafNodes, inputTree, *refPointIndexTree, *refSignFlagsTree,
5093 mQuantizedSeamPoints.get(), isovalue));
5098 const bool referenceMeshing = refSignFlagsTree && refPointIndexTree && refInputTree;
5103 boost::scoped_array<Index32> leafNodeOffsets(
new Index32[signFlagsLeafNodes.size()]);
5108 pointIndexLeafNodes, signFlagsLeafNodes, isovalue, adaptivityThreshold, invertSurfaceOrientation);
5110 if (mAdaptivityGrid && mAdaptivityGrid->type() == FloatGridType::gridType()) {
5111 const FloatGridType * adaptivityGrid =
static_cast<const FloatGridType*
>(mAdaptivityGrid.get());
5115 if (!adaptivityMask.empty()) {
5119 if (referenceMeshing) {
5123 tbb::parallel_for(auxiliaryLeafNodeRange, mergeOp);
5126 op(pointIndexLeafNodes, signFlagsLeafNodes, leafNodeOffsets);
5128 tbb::parallel_for(auxiliaryLeafNodeRange, op);
5133 op(signFlagsLeafNodes, leafNodeOffsets);
5135 tbb::parallel_for(auxiliaryLeafNodeRange, op);
5141 for (
size_t n = 0, N = signFlagsLeafNodes.size(); n < N; ++n) {
5142 const Index32 tmp = leafNodeOffsets[n];
5143 leafNodeOffsets[n] = pointCount;
5147 mPointListSize = size_t(pointCount);
5149 mPointFlags.clear();
5157 op(mPoints.get(), inputTree, pointIndexLeafNodes,
5158 signFlagsLeafNodes, leafNodeOffsets, transform, mIsovalue);
5160 if (referenceMeshing) {
5161 mPointFlags.resize(mPointListSize);
5162 op.setRefData(*refInputTree, *refPointIndexTree, *refSignFlagsTree,
5163 mQuantizedSeamPoints.get(), &mPointFlags.front());
5166 tbb::parallel_for(auxiliaryLeafNodeRange, op);
5172 mPolygonPoolListSize = signFlagsLeafNodes.size();
5173 mPolygons.reset(
new PolygonPool[mPolygonPoolListSize]);
5180 op(signFlagsLeafNodes, signFlagsTree, pointIndexTree, mPolygons, invertSurfaceOrientation);
5182 if (referenceMeshing) {
5186 tbb::parallel_for(auxiliaryLeafNodeRange, op);
5193 op(signFlagsLeafNodes, signFlagsTree, pointIndexTree, mPolygons, invertSurfaceOrientation);
5195 if (referenceMeshing) {
5199 tbb::parallel_for(auxiliaryLeafNodeRange, op);
5203 signFlagsTree.clear();
5204 pointIndexTree.clear();
5207 if (adaptive && mRelaxDisorientedTriangles) {
5209 inputTree, transform, mPolygons, mPolygonPoolListSize, mPoints, mPointListSize);
5213 if (referenceMeshing) {
5215 mPolygons, mPolygonPoolListSize, mPoints, mPointListSize, mPointFlags);
5227 template<
typename Gr
idType>
5228 inline typename boost::enable_if<boost::is_scalar<typename GridType::ValueType>,
void>::type
5230 const GridType& grid,
5231 std::vector<Vec3s>& points,
5232 std::vector<Vec3I>& triangles,
5233 std::vector<Vec4I>& quads,
5238 VolumeToMesh mesher(isovalue, adaptivity, relaxDisorientedTriangles);
5247 tbb::parallel_for(tbb::blocked_range<size_t>(0, points.size()), ptnCpy);
5254 size_t numQuads = 0, numTriangles = 0;
5256 openvdb::tools::PolygonPool& polygons = polygonPoolList[n];
5257 numTriangles += polygons.numTriangles();
5258 numQuads += polygons.numQuads();
5262 triangles.resize(numTriangles);
5264 quads.resize(numQuads);
5268 size_t qIdx = 0, tIdx = 0;
5270 openvdb::tools::PolygonPool& polygons = polygonPoolList[n];
5272 for (
size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
5273 quads[qIdx++] = polygons.quad(i);
5276 for (
size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
5277 triangles[tIdx++] = polygons.triangle(i);
5283 template<
typename Gr
idType>
5284 inline typename boost::disable_if<boost::is_scalar<typename GridType::ValueType>,
void>::type
5287 std::vector<Vec3s>&,
5288 std::vector<Vec3I>&,
5289 std::vector<Vec4I>&,
5298 template<
typename Gr
idType>
5301 const GridType& grid,
5302 std::vector<Vec3s>& points,
5303 std::vector<Vec3I>& triangles,
5304 std::vector<Vec4I>& quads,
5309 doVolumeToMesh(grid, points, triangles, quads, isovalue, adaptivity, relaxDisorientedTriangles);
5313 template<
typename Gr
idType>
5316 const GridType& grid,
5317 std::vector<Vec3s>& points,
5318 std::vector<Vec4I>& quads,
5321 std::vector<Vec3I> triangles;
5322 doVolumeToMesh(grid, points, triangles, quads, isovalue, 0.0,
true);
5333 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
Base class for typed trees.
Definition: Tree.h:65
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:374
OPENVDB_API const Index32 INVALID_IDX
Vec3< float > Vec3s
Definition: Vec3.h:650
boost::shared_ptr< const GridBase > ConstPtr
Definition: Grid.h:107
const Coord & origin() const
Return the grid index coordinates of this node's local origin.
Definition: LeafNode.h:475
T & z()
Definition: Vec3.h:99
Abstract base class for typed grids.
Definition: Grid.h:103
Vec4< int32_t > Vec4i
Definition: Vec4.h:581
boost::shared_ptr< const TreeBase > ConstPtr
Definition: Tree.h:69
Coord & reset(Int32 x, Int32 y, Int32 z)
Definition: Coord.h:94
static Coord offsetToLocalCoord(Index n)
Return the local coordinates for a linear table offset, where offset 0 has coordinates (0...
Definition: LeafNode.h:1286
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
boost::shared_ptr< TreeBase > Ptr
Definition: Tree.h:68
int getValueDepth(const Coord &xyz) const
Definition: ValueAccessor.h:275
Vec3< int32_t > Vec3i
Definition: Vec3.h:648
math::Transform & transform()
Return a reference to this grid's transform, which might be shared with other grids.
Definition: Grid.h:335
TreeType & tree()
Return a reference to this grid's tree, which might be shared with other grids.
Definition: Grid.h:761
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: ValueAccessor.h:256
Index32 Index
Definition: Types.h:58
T dot(const Vec3< T > &v) const
Dot product.
Definition: Vec3.h:203
int16_t Int16
Definition: Types.h:59
#define OPENVDB_VERSION_NAME
Definition: version.h:43
Vec3< double > Vec3d
Definition: Vec3.h:651
bool isValueOn(const Coord &xyz) const
Return the active state of the voxel at the given coordinates.
Definition: ValueAccessor.h:263
Mat3< double > Mat3d
Definition: Mat3.h:712
Templated block class to hold specific data types and a fixed number of values determined by Log2Dim...
Definition: LeafNode.h:65
Gradient operators defined in index space of various orders.
Definition: Operators.h:122
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:429
Definition: Exceptions.h:39
uint32_t Index32
Definition: Types.h:56
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: LeafNode.h:1311
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:803
ValueOnCIter cbeginValueOn() const
Definition: LeafNode.h:594
bool normalize(T eps=T(1.0e-7))
this = normalized this
Definition: Vec3.h:348
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:393
Definition: Exceptions.h:87
Definition: LeafNode.h:510
Definition: Exceptions.h:88
TreeType & tree() const
Return a reference to the tree associated with this accessor.
Definition: ValueAccessor.h:147
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
bool probeValue(const Coord &xyz, ValueType &value) const
Return the active state of the voxel as well as its value.
Definition: ValueAccessor.h:266
Mat3 transpose() const
returns transpose of this
Definition: Mat3.h:498
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:54
static const Index DIM
Definition: LeafNode.h:77
void setValue(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:287