31 #ifndef OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
32 #define OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
34 #include <openvdb/tree/ValueAccessor.h>
35 #include <openvdb/util/Util.h>
36 #include <openvdb/math/Operators.h>
37 #include <openvdb/tools/Morphology.h>
38 #include <openvdb/tree/LeafManager.h>
40 #include <boost/scoped_array.hpp>
41 #include <boost/scoped_ptr.hpp>
43 #include <tbb/blocked_range.h>
44 #include <tbb/parallel_for.h>
45 #include <tbb/parallel_reduce.h>
66 template<
typename Gr
idType>
70 std::vector<Vec3s>& points,
71 std::vector<Vec4I>& quads,
72 double isovalue = 0.0);
83 template<
typename Gr
idType>
87 std::vector<Vec3s>& points,
88 std::vector<Vec3I>& triangles,
89 std::vector<Vec4I>& quads,
90 double isovalue = 0.0,
91 double adaptivity = 0.0);
113 , mTriangleFlags(NULL)
121 mQuadFlags.reset(
new char[mNumQuads]);
128 mQuadFlags.reset(NULL);
133 mNumTriangles = size;
135 mTriangleFlags.reset(
new char[mNumTriangles]);
141 mTriangles.reset(NULL);
142 mTriangleFlags.reset(NULL);
145 const size_t&
numQuads()
const {
return mNumQuads; }
157 const char&
quadFlags(
size_t n)
const {
return mQuadFlags[n]; }
163 size_t mNumQuads, mNumTriangles;
164 boost::scoped_array<openvdb::Vec4I> mQuads;
165 boost::scoped_array<openvdb::Vec3I> mTriangles;
166 boost::scoped_array<char> mQuadFlags, mTriangleFlags;
186 VolumeToMesh(
double isovalue = 0,
double adaptivity = 0);
189 const size_t& pointListSize()
const;
193 const size_t& polygonPoolListSize()
const;
197 template<
typename Gr
idT>
198 void operator()(
const GridT&);
225 bool smoothSeams =
true);
247 void partition(
unsigned partitions = 1,
unsigned activePart = 0);
255 size_t mPointListSize, mPolygonPoolListSize;
256 double mIsovalue, mPrimAdaptivity, mSecAdaptivity;
261 TreeBase::Ptr mRefEdgeTree, mRefTopologyMaskTree, mSeamPointTree;
263 bool mSmoothSeams, mInvertSurfaceMask;
264 unsigned mPartitions, mActivePart;
278 const std::vector<Vec3d>& points,
279 const std::vector<Vec3d>& normals)
285 if (points.empty())
return avgPos;
287 for (
size_t n = 0, N = points.size(); n < N; ++n) {
291 avgPos /= double(points.size());
295 double m00=0,m01=0,m02=0,
302 for (
size_t n = 0, N = points.size(); n < N; ++n) {
304 const Vec3d& n_ref = normals[n];
307 m00 += n_ref[0] * n_ref[0];
308 m11 += n_ref[1] * n_ref[1];
309 m22 += n_ref[2] * n_ref[2];
311 m01 += n_ref[0] * n_ref[1];
312 m02 += n_ref[0] * n_ref[2];
313 m12 += n_ref[1] * n_ref[2];
316 rhs += n_ref * n_ref.
dot(points[n] - avgPos);
341 Mat3d D = Mat3d::identity();
345 tolerance =
std::max(tolerance, std::abs(eigenValues[2]));
349 for (
int i = 0; i < 3; ++i ) {
350 if (std::abs(eigenValues[i]) < tolerance) {
354 D[i][i] = 1.0 / eigenValues[i];
360 Mat3d pseudoInv = eigenVectors * D * eigenVectors.
transpose();
361 return avgPos + pseudoInv * rhs;
415 template<
class AccessorT>
417 typename AccessorT::ValueType isovalue,
int dim)
421 if (accessor.getValue(coord) < isovalue) signs |= 1u;
423 if (accessor.getValue(coord) < isovalue) signs |= 2u;
425 if (accessor.getValue(coord) < isovalue) signs |= 4u;
427 if (accessor.getValue(coord) < isovalue) signs |= 8u;
428 coord[1] += dim; coord[2] = ijk[2];
429 if (accessor.getValue(coord) < isovalue) signs |= 16u;
431 if (accessor.getValue(coord) < isovalue) signs |= 32u;
433 if (accessor.getValue(coord) < isovalue) signs |= 64u;
435 if (accessor.getValue(coord) < isovalue) signs |= 128u;
440 template<
class AccessorT>
443 typename AccessorT::ValueType isovalue,
const int dim)
449 p[0] = accessor.getValue(coord) < isovalue;
451 p[1] = accessor.getValue(coord) < isovalue;
453 p[2] = accessor.getValue(coord) < isovalue;
455 p[3] = accessor.getValue(coord) < isovalue;
456 coord[1] += dim; coord[2] = ijk[2];
457 p[4] = accessor.getValue(coord) < isovalue;
459 p[5] = accessor.getValue(coord) < isovalue;
461 p[6] = accessor.getValue(coord) < isovalue;
463 p[7] = accessor.getValue(coord) < isovalue;
467 if (p[0]) signs |= 1u;
468 if (p[1]) signs |= 2u;
469 if (p[2]) signs |= 4u;
470 if (p[3]) signs |= 8u;
471 if (p[4]) signs |= 16u;
472 if (p[5]) signs |= 32u;
473 if (p[6]) signs |= 64u;
474 if (p[7]) signs |= 128u;
480 int i = ijk[0], ip = ijk[0] + hDim, ipp = ijk[0] + dim;
481 int j = ijk[1], jp = ijk[1] + hDim, jpp = ijk[1] + dim;
482 int k = ijk[2], kp = ijk[2] + hDim, kpp = ijk[2] + dim;
485 coord.
reset(ip, j, k);
486 m = accessor.getValue(coord) < isovalue;
487 if (p[0] != m && p[1] != m)
return true;
490 coord.
reset(ipp, j, kp);
491 m = accessor.getValue(coord) < isovalue;
492 if (p[1] != m && p[2] != m)
return true;
495 coord.
reset(ip, j, kpp);
496 m = accessor.getValue(coord) < isovalue;
497 if (p[2] != m && p[3] != m)
return true;
500 coord.
reset(i, j, kp);
501 m = accessor.getValue(coord) < isovalue;
502 if (p[0] != m && p[3] != m)
return true;
505 coord.
reset(ip, jpp, k);
506 m = accessor.getValue(coord) < isovalue;
507 if (p[4] != m && p[5] != m)
return true;
510 coord.
reset(ipp, jpp, kp);
511 m = accessor.getValue(coord) < isovalue;
512 if (p[5] != m && p[6] != m)
return true;
515 coord.
reset(ip, jpp, kpp);
516 m = accessor.getValue(coord) < isovalue;
517 if (p[6] != m && p[7] != m)
return true;
520 coord.
reset(i, jpp, kp);
521 m = accessor.getValue(coord) < isovalue;
522 if (p[7] != m && p[4] != m)
return true;
525 coord.
reset(i, jp, k);
526 m = accessor.getValue(coord) < isovalue;
527 if (p[0] != m && p[4] != m)
return true;
530 coord.
reset(ipp, jp, k);
531 m = accessor.getValue(coord) < isovalue;
532 if (p[1] != m && p[5] != m)
return true;
535 coord.
reset(ipp, jp, kpp);
536 m = accessor.getValue(coord) < isovalue;
537 if (p[2] != m && p[6] != m)
return true;
541 coord.
reset(i, jp, kpp);
542 m = accessor.getValue(coord) < isovalue;
543 if (p[3] != m && p[7] != m)
return true;
549 coord.
reset(ip, jp, k);
550 m = accessor.getValue(coord) < isovalue;
551 if (p[0] != m && p[1] != m && p[4] != m && p[5] != m)
return true;
554 coord.
reset(ipp, jp, kp);
555 m = accessor.getValue(coord) < isovalue;
556 if (p[1] != m && p[2] != m && p[5] != m && p[6] != m)
return true;
559 coord.
reset(ip, jp, kpp);
560 m = accessor.getValue(coord) < isovalue;
561 if (p[2] != m && p[3] != m && p[6] != m && p[7] != m)
return true;
564 coord.
reset(i, jp, kp);
565 m = accessor.getValue(coord) < isovalue;
566 if (p[0] != m && p[3] != m && p[4] != m && p[7] != m)
return true;
569 coord.
reset(ip, j, kp);
570 m = accessor.getValue(coord) < isovalue;
571 if (p[0] != m && p[1] != m && p[2] != m && p[3] != m)
return true;
574 coord.
reset(ip, jpp, kp);
575 m = accessor.getValue(coord) < isovalue;
576 if (p[4] != m && p[5] != m && p[6] != m && p[7] != m)
return true;
579 coord.
reset(ip, jp, kp);
580 m = accessor.getValue(coord) < isovalue;
581 if (p[0] != m && p[1] != m && p[2] != m && p[3] != m &&
582 p[4] != m && p[5] != m && p[6] != m && p[7] != m)
return true;
591 template <
class LeafType>
595 Coord ijk, end = start;
600 for (ijk[0] = start[0]; ijk[0] < end[0]; ++ijk[0]) {
601 for (ijk[1] = start[1]; ijk[1] < end[1]; ++ijk[1]) {
602 for (ijk[2] = start[2]; ijk[2] < end[2]; ++ijk[2]) {
603 if(leaf.isValueOn(ijk)) leaf.setValue(ijk, regionId);
612 template <
class LeafType>
615 typename LeafType::ValueType::value_type adaptivity)
617 if (adaptivity < 1e-6)
return false;
619 typedef typename LeafType::ValueType VecT;
620 Coord ijk, end = start;
625 std::vector<VecT> norms;
626 for (ijk[0] = start[0]; ijk[0] < end[0]; ++ijk[0]) {
627 for (ijk[1] = start[1]; ijk[1] < end[1]; ++ijk[1]) {
628 for (ijk[2] = start[2]; ijk[2] < end[2]; ++ijk[2]) {
630 if(!leaf.isValueOn(ijk))
continue;
631 norms.push_back(leaf.getValue(ijk));
636 size_t N = norms.size();
637 for (
size_t ni = 0; ni < N; ++ni) {
638 VecT n_i = norms[ni];
639 for (
size_t nj = 0; nj < N; ++nj) {
640 VecT n_j = norms[nj];
641 if ((1.0 - n_i.dot(n_j)) > adaptivity)
return false;
651 template <
class TreeT>
655 typedef std::vector<const typename TreeT::LeafNodeType *>
ListT;
659 mLeafNodes.reserve(tree.leafCount());
660 typename TreeT::LeafCIter iter = tree.cbeginLeaf();
661 for ( ; iter; ++iter) mLeafNodes.push_back(iter.getLeaf());
664 size_t size()
const {
return mLeafNodes.size(); }
666 const typename TreeT::LeafNodeType*
operator[](
size_t n)
const
667 {
return mLeafNodes[n]; }
670 {
return tbb::blocked_range<size_t>(0, mLeafNodes.size()); }
679 template <
class TreeT>
683 typedef std::vector<typename TreeT::LeafNodeType *>
ListT;
687 mLeafNodes.reserve(tree.leafCount());
688 typename TreeT::LeafIter iter = tree.beginLeaf();
689 for ( ; iter; ++iter) mLeafNodes.push_back(iter.getLeaf());
692 size_t size()
const {
return mLeafNodes.size(); }
695 {
return mLeafNodes[n]; }
698 {
return tbb::blocked_range<size_t>(0, mLeafNodes.size()); }
709 template<
typename DistTreeT>
713 typedef typename DistTreeT::template ValueConverter<char>::Type
CharTreeT;
714 typedef typename DistTreeT::template ValueConverter<bool>::Type
BoolTreeT;
715 typedef typename DistTreeT::template ValueConverter<Vec3s>::Type
Vec3sTreeT;
720 , mTopologyMaskTree(NULL)
721 , mSeamPointTree(NULL)
722 , mSeamMaskTree(typename
BoolTreeT::Ptr())
723 , mSmoothingMaskTree(typename
BoolTreeT::Ptr())
730 return mDistTree && mEdgeTree && mTopologyMaskTree && mSeamMaskTree;
745 template <
class DistTreeT>
755 inline void operator()(
const tbb::blocked_range<size_t>&)
const;
758 std::vector<size_t>& mLeafRegionCount;
762 template <
class DistTreeT>
765 std::vector<size_t>& leafRegionCount)
767 , mLeafRegionCount(leafRegionCount)
772 template <
class DistTreeT>
775 : mLeafNodes(rhs.mLeafNodes)
776 , mLeafRegionCount(rhs.mLeafRegionCount)
781 template <
class DistTreeT>
785 tbb::parallel_for(mLeafNodes.getRange(), *
this);
789 template <
class DistTreeT>
793 (*this)(mLeafNodes.getRange());
797 template <
class DistTreeT>
801 for (
size_t n = range.begin(); n != range.end(); ++n) {
802 mLeafRegionCount[n] = size_t(mLeafNodes[n]->onVoxelCount());
810 template <
class DistTreeT>
815 typedef typename DistTreeT::template ValueConverter<bool>::Type
BoolTreeT;
816 typedef typename DistTreeT::template ValueConverter<int>::Type
IntTreeT;
817 typedef typename DistTreeT::template ValueConverter<float>::Type
FloatTreeT;
821 const DistTreeT& distTree,
823 std::vector<size_t>& leafRegionCount,
835 void run(
bool threaded =
true);
837 void operator()(
const tbb::blocked_range<size_t>&)
const;
840 const DistTreeT& mDistTree;
842 std::vector<size_t>& mLeafRegionCount;
852 template <
class DistTreeT>
854 const DistTreeT& distTree,
856 std::vector<size_t>& leafRegionCount,
860 : mDistTree(distTree)
861 , mAuxLeafs(auxLeafs)
862 , mLeafRegionCount(leafRegionCount)
864 , mAdaptivity(adaptivity)
867 , mAdaptivityGrid(NULL)
868 , mPointMask(pointMask)
873 template <
class DistTreeT>
876 : mDistTree(rhs.mDistTree)
877 , mAuxLeafs(rhs.mAuxLeafs)
878 , mLeafRegionCount(rhs.mLeafRegionCount)
879 , mIsovalue(rhs.mIsovalue)
880 , mAdaptivity(rhs.mAdaptivity)
881 , mRefData(rhs.mRefData)
882 , mTransform(rhs.mTransform)
883 , mAdaptivityGrid(rhs.mAdaptivityGrid)
884 , mPointMask(rhs.mPointMask)
889 template <
class DistTreeT>
894 tbb::parallel_for(mAuxLeafs.getRange(), *
this);
896 (*this)(mAuxLeafs.getRange());
901 template <
class DistTreeT>
909 template <
class DistTreeT>
914 mTransform = &distGridXForm;
915 mAdaptivityGrid = &adaptivityField;
919 template <
class DistTreeT>
924 typedef typename BoolTreeT::LeafNodeType BoolLeafT;
925 typedef typename IntTreeT::LeafNodeType IntLeafT;
926 typedef typename BoolLeafT::template ValueConverter<Vec3T>::Type Vec3LeafT;
927 typedef typename BoolLeafT::template ValueConverter<DistValueT>::Type DistLeafT;
929 typedef typename IntLeafT::ValueOnIter IntIterT;
930 typedef typename BoolLeafT::ValueOnCIter BoolCIterT;
936 boost::scoped_ptr<FloatTreeCAccessorT> adaptivityAcc;
937 if (mAdaptivityGrid) {
938 adaptivityAcc.reset(
new FloatTreeCAccessorT(mAdaptivityGrid->tree()));
941 boost::scoped_ptr<BoolTreeAccessorT> seamMaskAcc;
942 boost::scoped_ptr<BoolTreeCAccessorT> topologyMaskAcc;
943 if (mRefData && mRefData->isValid()) {
944 seamMaskAcc.reset(
new BoolTreeAccessorT(*mRefData->mSeamMaskTree.get()));
945 topologyMaskAcc.reset(
new BoolTreeCAccessorT(*mRefData->mTopologyMaskTree));
947 const bool hasRefData = seamMaskAcc && topologyMaskAcc;
949 const int LeafDim = BoolLeafT::DIM;
954 Vec3LeafT gradientBuffer;
955 Coord ijk, coord, end;
957 for (
size_t n = range.begin(); n != range.end(); ++n) {
960 IntLeafT& auxLeaf = *mAuxLeafs[n];
962 const Coord& origin = auxLeaf.getOrigin();
963 end[0] = origin[0] + LeafDim;
964 end[1] = origin[1] + LeafDim;
965 end[2] = origin[2] + LeafDim;
971 const BoolLeafT* seamMask = seamMaskAcc->probeConstLeaf(origin);
972 if (seamMask != NULL) {
973 for (BoolCIterT it = seamMask->cbeginValueOn(); it; ++it) {
975 coord[0] = ijk[0] - (ijk[0] % 2);
976 coord[1] = ijk[1] - (ijk[1] % 2);
977 coord[2] = ijk[2] - (ijk[2] % 2);
978 mask.setActiveState(coord,
true);
981 if (topologyMaskAcc->probeConstLeaf(origin) == NULL) {
982 adaptivity = mRefData->mInternalAdaptivity;
987 DistLeafT adaptivityLeaf(origin, adaptivity);
989 if (mAdaptivityGrid) {
990 for (
Index offset = 0; offset < DistLeafT::NUM_VALUES; ++offset) {
991 ijk = adaptivityLeaf.offset2globalCoord(offset);
992 Vec3d xyz = mAdaptivityGrid->transform().worldToIndex(
993 mTransform->indexToWorld(ijk));
995 adaptivityLeaf.setValueOnly(offset, tmpA * adaptivity);
1000 for (IntIterT it = auxLeaf.beginValueOn(); it; ++it) {
1001 ijk = it.getCoord();
1002 coord[0] = ijk[0] - (ijk[0] % 2);
1003 coord[1] = ijk[1] - (ijk[1] % 2);
1004 coord[2] = ijk[2] - (ijk[2] % 2);
1005 if(mask.isValueOn(coord))
continue;
1006 mask.setActiveState(coord,
isAmbiguous(distAcc, ijk, mIsovalue, 1));
1011 for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
1012 for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
1013 for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
1015 mask.setActiveState(ijk,
true);
1022 gradientBuffer.setValuesOff();
1024 for (IntIterT it = auxLeaf.beginValueOn(); it; ++it) {
1026 ijk = it.getCoord();
1027 coord[0] = ijk[0] - (ijk[0] % dim);
1028 coord[1] = ijk[1] - (ijk[1] % dim);
1029 coord[2] = ijk[2] - (ijk[2] % dim);
1030 if(mask.isValueOn(coord))
continue;
1038 gradientBuffer.setValue(ijk, norm);
1041 int regionId = 1, next_dim = dim << 1;
1044 for (ijk[0] = 0; ijk[0] < LeafDim; ijk[0] += dim) {
1045 coord[0] = ijk[0] - (ijk[0] % next_dim);
1046 for (ijk[1] = 0; ijk[1] < LeafDim; ijk[1] += dim) {
1047 coord[1] = ijk[1] - (ijk[1] % next_dim);
1048 for (ijk[2] = 0; ijk[2] < LeafDim; ijk[2] += dim) {
1049 coord[2] = ijk[2] - (ijk[2] % next_dim);
1050 adaptivity = adaptivityLeaf.getValue(ijk);
1051 if(mask.isValueOn(ijk) || !
isMergable(gradientBuffer, ijk, dim, adaptivity)) {
1052 mask.setActiveState(coord,
true);
1062 for (dim = 4; dim < LeafDim; dim = dim << 1) {
1063 next_dim = dim << 1;
1064 coord[0] = ijk[0] - (ijk[0] % next_dim);
1065 for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
1066 coord[1] = ijk[1] - (ijk[1] % next_dim);
1067 for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
1068 coord[2] = ijk[2] - (ijk[2] % next_dim);
1069 for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
1070 adaptivity = adaptivityLeaf.getValue(ijk);
1071 if (mask.isValueOn(ijk) ||
isNonManifold(distAcc, ijk, mIsovalue, dim) ||
1072 !
isMergable(gradientBuffer, ijk, dim, adaptivity)) {
1073 mask.setActiveState(coord,
true);
1082 adaptivity = adaptivityLeaf.getValue(origin);
1083 if (!(mask.isValueOn(origin) ||
isNonManifold(distAcc, origin, mIsovalue, LeafDim))
1084 &&
isMergable(gradientBuffer, origin, LeafDim, adaptivity)) {
1085 mergeVoxels(auxLeaf, origin, LeafDim, regionId++);
1089 const BoolLeafT * ptnMaskLeaf = mPointMask->probeConstLeaf(auxLeaf.getOrigin());
1093 for (IntIterT it = auxLeaf.beginValueOn(); it; ++it) {
1094 if (it.getValue() == 0 && !ptnMaskLeaf->isValueOn(it.pos())) {
1099 IntLeafT tmpLeaf(auxLeaf);
1101 while (tmpLeaf.onVoxelCount() > 0) {
1103 IntIterT it = tmpLeaf.beginValueOn();
1104 regionId = it.getValue();
1105 bool removeRegion =
true;
1107 if (it.getValue() == regionId) {
1108 if (ptnMaskLeaf->isValueOn(it.pos())) {
1109 removeRegion =
false;
1116 it = auxLeaf.beginValueOn();
1118 if (it.getValue() == regionId) {
1125 mLeafRegionCount[n] = 0;
1131 size_t numVoxels = 0;
1132 IntLeafT tmpLeaf(auxLeaf);
1133 for (IntIterT it = tmpLeaf.beginValueOn(); it; ++it) {
1134 if(it.getValue() == 0) {
1140 while (tmpLeaf.onVoxelCount() > 0) {
1142 IntIterT it = tmpLeaf.beginValueOn();
1143 regionId = it.getValue();
1145 if (it.getValue() == regionId) it.setValueOff();
1149 mLeafRegionCount[n] = numVoxels;
1157 template <
class DistTreeT>
1163 typedef typename DistTreeT::template ValueConverter<int>::Type
IntTreeT;
1166 const DistTreeT& distTree,
1168 std::vector<size_t>& leafIndices,
1169 const openvdb::math::Transform& xform,
1180 void operator()(
const tbb::blocked_range<size_t>&)
const;
1183 const DistTreeT& mDistTree;
1185 const std::vector<size_t>& mLeafIndices;
1186 const openvdb::math::Transform& mTransform;
1188 const double mIsovalue;
1192 double root(
double v0,
double v1)
const {
return (mIsovalue - v0) / (v1 - v0); }
1197 template <
class DistTreeT>
1199 const DistTreeT& distTree,
1201 std::vector<size_t>& leafIndices,
1202 const openvdb::math::Transform& xform,
1205 : mDistTree(distTree)
1206 , mAuxLeafs(auxLeafs)
1207 , mLeafIndices(leafIndices)
1216 template <
class DistTreeT>
1218 : mDistTree(rhs.mDistTree)
1219 , mAuxLeafs(rhs.mAuxLeafs)
1220 , mLeafIndices(rhs.mLeafIndices)
1221 , mTransform(rhs.mTransform)
1222 , mPoints(rhs.mPoints)
1223 , mIsovalue(rhs.mIsovalue)
1224 , mRefData(rhs.mRefData)
1229 template <
class DistTreeT>
1234 mRefData = &refData;
1237 template <
class DistTreeT>
1241 tbb::parallel_for(mAuxLeafs.getRange(), *
this);
1245 template <
class DistTreeT>
1249 (*this)(mAuxLeafs.getRange());
1253 template <
class DistTreeT>
1257 typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
1258 typedef typename DistTreeT::template ValueConverter<Vec3s>::Type Vec3sTreeT;
1263 typedef typename BoolTreeT::LeafNodeType BoolLeafT;
1264 typedef typename IntTreeT::LeafNodeType IntLeafT;
1265 typedef typename Vec3sTreeT::LeafNodeType Vec3sLeafT;
1267 boost::scoped_ptr<DistTreeAccessorT> refDistAcc;
1268 boost::scoped_ptr<BoolTreeAccessorT> refMaskAcc, refSmoothMaskAcc;
1269 boost::scoped_ptr<Vec3sTreeAccessorT> refPtnAcc;
1271 if (mRefData && mRefData->isValid()) {
1273 refMaskAcc.reset(
new BoolTreeAccessorT(*mRefData->mTopologyMaskTree));
1274 refSmoothMaskAcc.reset(
new BoolTreeAccessorT(*mRefData->mSmoothingMaskTree));
1275 refPtnAcc.reset(
new Vec3sTreeAccessorT(*mRefData->mSeamPointTree));
1279 const bool hasRefData = refDistAcc && refMaskAcc;
1280 typename IntTreeT::LeafNodeType::ValueOnIter auxIter;
1286 for (
size_t n = range.begin(); n != range.end(); ++n) {
1288 size_t idx = mLeafIndices[n];
1289 IntLeafT& auxLeaf = *mAuxLeafs[n];
1291 BoolLeafT* maskLeaf = NULL;
1292 BoolLeafT* smoothMaskLeaf = NULL;
1293 Vec3sLeafT* ptnLeaf = NULL;
1295 maskLeaf = refMaskAcc->probeLeaf(auxLeaf.getOrigin());
1296 smoothMaskLeaf = refSmoothMaskAcc->probeLeaf(auxLeaf.getOrigin());
1297 ptnLeaf = refPtnAcc->probeLeaf(auxLeaf.getOrigin());
1300 for (auxIter = auxLeaf.beginValueOn(); auxIter; ++auxIter) {
1302 if(auxIter.getValue() == 0) {
1304 auxIter.setValue(idx);
1305 auxIter.setValueOff();
1306 ijk = auxIter.getCoord();
1308 if (hasRefData && maskLeaf && maskLeaf->isValueOn(ijk)) {
1310 if (ptnLeaf && ptnLeaf->isValueOn(ijk)) {
1311 avg = ptnLeaf->getValue(ijk);
1313 int e1 = calcAvgPoint(*refDistAcc.get(), ijk, avg);
1316 int e2 = calcAvgPoint(distAcc, ijk, tmp);
1317 if((e2 & (~e1)) != 0) smoothMaskLeaf->setValueOn(ijk);
1321 calcAvgPoint(distAcc, ijk, avg);
1325 ptn[0] = float(avg[0]);
1326 ptn[1] = float(avg[1]);
1327 ptn[2] = float(avg[2]);
1333 while(auxLeaf.onVoxelCount() > 0) {
1339 auxIter = auxLeaf.beginValueOn();
1340 int regionId = auxIter.getValue(), points = 0;
1342 for (; auxIter; ++auxIter) {
1343 if(auxIter.getValue() == regionId) {
1345 auxIter.setValue(idx);
1346 auxIter.setValueOff();
1347 ijk = auxIter.getCoord();
1349 if (hasRefData && maskLeaf && maskLeaf->isValueOn(ijk)) {
1350 calcAvgPoint(*refDistAcc.get(), ijk, tmp);
1352 calcAvgPoint(distAcc, ijk, tmp);
1362 double w = 1.0 / double(points);
1369 ptn[0] = float(avg[0]);
1370 ptn[1] = float(avg[1]);
1371 ptn[2] = float(avg[2]);
1377 template <
class DistTreeT>
1388 values[0] = double(acc.getValue(coord));
1391 values[1] = double(acc.getValue(coord));
1394 values[2] = double(acc.getValue(coord));
1397 values[3] = double(acc.getValue(coord));
1399 coord[1] += 1; coord[2] = ijk[2];
1400 values[4] = double(acc.getValue(coord));
1403 values[5] = double(acc.getValue(coord));
1406 values[6] = double(acc.getValue(coord));
1409 values[7] = double(acc.getValue(coord));
1412 for (
int n = 0; n < 8; ++n) signMask[n] = (values[n] < mIsovalue);
1414 int samples = 0, edgeFlags = 0;
1419 if (signMask[0] != signMask[1]) {
1420 avg[0] += root(values[0], values[1]);
1425 if (signMask[1] != signMask[2]) {
1427 avg[2] += root(values[1], values[2]);
1432 if (signMask[3] != signMask[2]) {
1433 avg[0] += root(values[3], values[2]);
1439 if (signMask[0] != signMask[3]) {
1440 avg[2] += root(values[0], values[3]);
1445 if (signMask[4] != signMask[5]) {
1446 avg[0] += root(values[4], values[5]);
1452 if (signMask[5] != signMask[6]) {
1455 avg[2] += root(values[5], values[6]);
1460 if (signMask[7] != signMask[6]) {
1461 avg[0] += root(values[7], values[6]);
1468 if (signMask[4] != signMask[7]) {
1470 avg[2] += root(values[4], values[7]);
1475 if (signMask[0] != signMask[4]) {
1476 avg[1] += root(values[0], values[4]);
1481 if (signMask[1] != signMask[5]) {
1483 avg[1] += root(values[1], values[5]);
1488 if (signMask[2] != signMask[6]) {
1490 avg[1] += root(values[2], values[6]);
1496 if (signMask[3] != signMask[7]) {
1497 avg[1] += root(values[3], values[7]);
1504 double w = 1.0 / double(samples);
1511 avg[0] += double(ijk[0]);
1512 avg[1] += double(ijk[1]);
1513 avg[2] += double(ijk[2]);
1515 avg = mTransform.indexToWorld(avg);
1530 mPolygonPool = &quadPool;
1538 mPolygonPool->
quad(mIdx) = verts;
1560 AdaptiveMeshOp(): mQuadIdx(0), mTriangleIdx(0), mPolygonPool(NULL), mTmpPolygonPool() {}
1564 mPolygonPool = &polygonPool;
1575 if (verts[0] != verts[1] && verts[0] != verts[2] && verts[0] != verts[3]
1576 && verts[1] != verts[2] && verts[1] != verts[3] && verts[2] != verts[3]) {
1577 mTmpPolygonPool.
quadFlags(mQuadIdx) = flags;
1578 addQuad(verts, reverse);
1580 verts[0] == verts[3] &&
1581 verts[1] != verts[2] &&
1582 verts[1] != verts[0] &&
1583 verts[2] != verts[0]) {
1585 addTriangle(verts[0], verts[1], verts[2], reverse);
1587 verts[1] == verts[2] &&
1588 verts[0] != verts[3] &&
1589 verts[0] != verts[1] &&
1590 verts[3] != verts[1]) {
1592 addTriangle(verts[0], verts[1], verts[3], reverse);
1594 verts[0] == verts[1] &&
1595 verts[2] != verts[3] &&
1596 verts[2] != verts[0] &&
1597 verts[3] != verts[0]) {
1599 addTriangle(verts[0], verts[2], verts[3], reverse);
1601 verts[2] == verts[3] &&
1602 verts[0] != verts[1] &&
1603 verts[0] != verts[2] &&
1604 verts[1] != verts[2]) {
1606 addTriangle(verts[0], verts[1], verts[2], reverse);
1614 for (
size_t i = 0; i < mQuadIdx; ++i) {
1615 mPolygonPool->
quad(i) = mTmpPolygonPool.
quad(i);
1621 for (
size_t i = 0; i < mTriangleIdx; ++i) {
1630 void addQuad(
const Vec4I& verts,
bool reverse)
1633 mTmpPolygonPool.
quad(mQuadIdx) = verts;
1635 Vec4I& quad = mTmpPolygonPool.
quad(mQuadIdx);
1644 void addTriangle(
unsigned v0,
unsigned v1,
unsigned v2,
bool reverse)
1660 size_t mQuadIdx, mTriangleIdx;
1661 PolygonPool *mPolygonPool;
1662 PolygonPool mTmpPolygonPool;
1669 template<
class DistTreeT,
class MeshingOp>
1673 typedef typename DistTreeT::template ValueConverter<char>::Type
CharTreeT;
1674 typedef typename DistTreeT::template ValueConverter<int>::Type
IntTreeT;
1684 void operator()(
const tbb::blocked_range<size_t>&)
const;
1695 template<
class DistTreeT,
class MeshingOp>
1698 : mEdgeLeafs(edgeLeafs)
1700 , mPolygonPoolList(polygonPoolList)
1706 template<
class DistTreeT,
class MeshingOp>
1708 : mEdgeLeafs(rhs.mEdgeLeafs)
1709 , mAuxTree(rhs.mAuxTree)
1710 , mPolygonPoolList(rhs.mPolygonPoolList)
1711 , mRefData(rhs.mRefData)
1716 template<
class DistTreeT,
class MeshingOp>
1721 mRefData = &refData;
1725 template<
class DistTreeT,
class MeshingOp>
1729 tbb::parallel_for(mEdgeLeafs.getRange(), *
this);
1733 template<
class DistTreeT,
class MeshingOp>
1737 (*this)(mEdgeLeafs.getRange());
1741 template<
class DistTreeT,
class MeshingOp>
1744 const tbb::blocked_range<size_t>& range)
const
1746 typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
1747 typedef typename BoolTreeT::LeafNodeType BoolLeafT;
1748 typedef typename CharTreeT::LeafNodeType CharLeafT;
1750 typedef openvdb::tree::ValueAccessor<const CharTreeT> CharTreeAccessorT;
1751 typedef openvdb::tree::ValueAccessor<const IntTreeT> IntTreeAccessorT;
1752 typedef openvdb::tree::ValueAccessor<const BoolTreeT> BoolTreeAccessorT;
1754 boost::scoped_ptr<CharTreeAccessorT> refEdgeAcc;
1755 boost::scoped_ptr<BoolTreeAccessorT> refMaskAcc;
1756 if (mRefData && mRefData->isValid()) {
1757 refEdgeAcc.reset(
new CharTreeAccessorT(*mRefData->mEdgeTree));
1758 refMaskAcc.reset(
new BoolTreeAccessorT(*mRefData->mSeamMaskTree.get()));
1760 const bool hasRefData = refEdgeAcc && refMaskAcc;
1763 typename CharTreeT::LeafNodeType::ValueOnCIter iter;
1764 IntTreeAccessorT auxAcc(mAuxTree);
1767 char refEdgeFlags, isSemLinePoly;
1774 for (
size_t n = range.begin(); n != range.end(); ++n) {
1776 const Coord origin = mEdgeLeafs[n]->getOrigin();
1780 iter = mEdgeLeafs[n]->cbeginValueOn();
1781 for (; iter; ++iter) {
1782 char edgeFlags = iter.getValue() >> 1;
1783 edgeCount += edgeFlags & 0x1;
1785 edgeFlags = edgeFlags >> 1;
1786 edgeCount += edgeFlags & 0x1;
1788 edgeFlags = edgeFlags >> 1;
1789 edgeCount += edgeFlags & 0x1;
1792 mesher.init(edgeCount, mPolygonPoolList[n]);
1794 const CharLeafT* refEdgeLeaf = NULL;
1795 const BoolLeafT* refMaskLeaf = NULL;
1798 refEdgeLeaf = refEdgeAcc->probeConstLeaf(origin);
1799 refMaskLeaf = refMaskAcc->probeConstLeaf(origin);
1802 iter = mEdgeLeafs[n]->cbeginValueOn();
1803 for (; iter; ++iter) {
1804 ijk = iter.getCoord();
1805 const char& edgeFlags = iter.getValue();
1807 const bool isInside = edgeFlags &
INSIDE;
1812 if(refEdgeLeaf) refEdgeFlags = refEdgeLeaf->getValue(ijk);
1813 if (refMaskLeaf && refMaskLeaf->isValueOn(ijk)) {
1819 int v0 = auxAcc.getValue(ijk);
1821 if (edgeFlags &
XEDGE) {
1824 coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2];
1825 quad[1] = auxAcc.getValue(coord);
1827 quad[2] = auxAcc.getValue(coord);
1829 quad[3] = auxAcc.getValue(coord);
1831 mesher.addPrim(quad, isInside,
1832 (isSemLinePoly | isExteriorPoly[
bool(refEdgeFlags & XEDGE)]));
1836 if (edgeFlags &
YEDGE) {
1839 coord[0] = ijk[0]; coord[1] = ijk[1]; coord[2] = ijk[2]-1;
1840 quad[1] = auxAcc.getValue(coord);
1842 quad[2] = auxAcc.getValue(coord);
1844 quad[3] = auxAcc.getValue(coord);
1846 mesher.addPrim(quad, isInside,
1847 (isSemLinePoly | isExteriorPoly[
bool(refEdgeFlags & YEDGE)]));
1850 if (edgeFlags &
ZEDGE) {
1853 coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2];
1854 quad[1] = auxAcc.getValue(coord);
1856 quad[2] = auxAcc.getValue(coord);
1858 quad[3] = auxAcc.getValue(coord);
1860 mesher.addPrim(quad, !isInside,
1861 (isSemLinePoly | isExteriorPoly[
bool(refEdgeFlags & ZEDGE)]));
1877 PartOp(
size_t leafCount,
size_t partitions,
size_t activePart)
1879 size_t leafSegments = leafCount / partitions;
1880 mStart = leafSegments * activePart;
1881 mEnd = activePart >= (partitions - 1) ? leafCount : mStart + leafSegments;
1884 template <
typename LeafNodeType>
1887 if (leafIndex < mStart || leafIndex >= mEnd) leaf.setValuesOff();
1891 size_t mStart, mEnd;
1894 template<
typename SrcTreeT>
1899 typedef typename SrcTreeT::template ValueConverter<bool>::Type
BoolTreeT;
1907 void run(
bool threaded =
true);
1915 void operator()(
const tbb::blocked_range<size_t>&);
1921 size_t mStart, mEnd;
1924 template<
typename SrcTreeT>
1926 : mLeafManager(leafs)
1930 size_t leafSegments = leafCount / partitions;
1931 mStart = leafSegments * activePart;
1932 mEnd = activePart >= (partitions - 1) ? leafCount : mStart + leafSegments;
1935 template<
typename SrcTreeT>
1937 : mLeafManager(rhs.mLeafManager)
1939 , mStart(rhs.mStart)
1945 template<
typename SrcTreeT>
1950 tbb::parallel_reduce(mLeafManager.getRange(), *
this);
1952 (*this)(mLeafManager.getRange());
1957 template<
typename SrcTreeT>
1964 typedef typename BoolTreeT::LeafNodeType BoolLeafT;
1965 typename SrcTreeT::LeafNodeType::ValueOnCIter iter;
1967 for (
size_t n = range.begin(); n != range.end(); ++n) {
1968 if (n < mStart || n >= mEnd)
continue;
1969 BoolLeafT* leaf = acc.
touchLeaf(mLeafManager.leaf(n).getOrigin());
1970 leaf->topologyUnion(mLeafManager.leaf(n));
1975 template<
typename SrcTreeT>
1980 typedef typename SrcTreeT::template ValueConverter<bool>::Type
BoolTreeT;
1992 void run(
bool threaded =
true);
2001 void operator()(
const tbb::blocked_range<size_t>&);
2014 template<
typename SrcTreeT>
2018 , mLeafManager(srcLeafs)
2019 , mSrcXForm(srcXForm)
2020 , mInvertMask(invertMask)
2025 template<
typename SrcTreeT>
2028 , mLeafManager(rhs.mLeafManager)
2029 , mSrcXForm(rhs.mSrcXForm)
2030 , mInvertMask(rhs.mInvertMask)
2036 template<
typename SrcTreeT>
2041 tbb::parallel_reduce(mLeafManager.getRange(), *
this);
2043 (*this)(mLeafManager.getRange());
2048 template<
typename SrcTreeT>
2054 typedef typename BoolTreeT::LeafNodeType BoolLeafT;
2059 typename SrcTreeT::LeafNodeType::ValueOnCIter iter;
2060 for (
size_t n = range.begin(); n != range.end(); ++n) {
2062 ijk = mLeafManager.leaf(n).getOrigin();
2063 BoolLeafT* leaf =
new BoolLeafT(ijk,
false);
2064 bool addLeaf =
false;
2066 if (maskXForm == mSrcXForm) {
2068 leaf->topologyUnion(mLeafManager.leaf(n));
2069 const BoolLeafT* maskLeaf = maskAcc.probeConstLeaf(ijk);
2071 leaf->topologyDifference(*maskLeaf,
false);
2072 addLeaf = !leaf->isEmpty();
2074 addLeaf = maskAcc.isValueOn(ijk) != mInvertMask;
2078 for (iter = mLeafManager.leaf(n).cbeginValueOn(); iter; ++iter) {
2079 ijk = iter.getCoord();
2080 xyz = maskXForm.
worldToIndex(mSrcXForm.indexToWorld(ijk));
2082 leaf->setValueOn(iter.pos());
2088 if (addLeaf) acc.
addLeaf(leaf);
2093 template<
typename IntTreeT>
2100 : mAuxAcc(aux), mBoundaryAuxAcc(boundaryAux) {}
2102 template <
typename LeafNodeType>
2105 Coord ijk = leaf.getOrigin();
2110 if (auxLeaf && boundaryLeaf) auxLeaf->merge(*boundaryLeaf);
2117 template<
typename SrcTreeT>
2121 typedef typename SrcTreeT::template ValueConverter<int>::Type
IntTreeT;
2122 typedef typename SrcTreeT::template ValueConverter<bool>::Type
BoolTreeT;
2133 void run(
bool threaded =
true);
2141 void operator()(
const tbb::blocked_range<size_t>&);
2148 const SrcTreeT& mSrcTree;
2153 template<
typename SrcTreeT>
2156 : mLeafManager(auxLeafs)
2158 , mMaskTree(boolTree)
2163 template<
typename SrcTreeT>
2165 : mLeafManager(rhs.mLeafManager)
2166 , mSrcTree(rhs.mSrcTree)
2167 , mMaskTree(rhs.mMaskTree)
2172 template<
typename SrcTreeT>
2177 tbb::parallel_reduce(mLeafManager.getRange(), *
this);
2179 (*this)(mLeafManager.getRange());
2184 template<
typename SrcTreeT>
2189 typedef typename BoolTreeT::LeafNodeType BoolLeafT;
2190 typedef typename SrcTreeT::LeafNodeType SrcLeafT;
2196 typename IntTreeT::LeafNodeType::ValueOnCIter iter;
2197 for (
size_t n = range.begin(); n != range.end(); ++n) {
2199 ijk = mLeafManager.leaf(n).getOrigin();
2202 if (!srcLeaf)
continue;
2204 BoolLeafT* leaf =
new BoolLeafT(ijk,
false);
2205 leaf->topologyUnion(*srcLeaf);
2210 leaf->topologyDifference(*maskLeaf,
false);
2212 if (!leaf->isEmpty()) acc.
addLeaf(leaf);
2228 typename AuxDataT =
int>
2233 typedef typename SrcTreeT::template ValueConverter<char>::Type
CharTreeT;
2234 typedef typename SrcTreeT::template ValueConverter<AuxDataT>::Type
AuxTreeT;
2242 AuxiliaryData(
const SrcTreeT&,
const LeafManagerT&,
double iso,
bool extraCheck =
false);
2244 void run(
bool threaded =
true);
2246 typename CharTreeT::Ptr
edgeTree()
const {
return mEdgeTree; }
2247 typename AuxTreeT::Ptr
auxTree()
const {
return mAuxTree; }
2254 void operator()(
const tbb::blocked_range<size_t>&);
2258 mEdgeTree->merge(*rhs.mEdgeTree);
2259 mAuxTree->merge(*rhs.mAuxTree);
2264 int edgeCheck(
const Coord& ijk,
const bool thisInside);
2266 const LeafManagerT& mLeafManager;
2268 const SrcTreeT& mSourceTree;
2271 typename CharTreeT::Ptr mEdgeTree;
2274 typename AuxTreeT::Ptr mAuxTree;
2277 const double mIsovalue;
2278 const bool mExtraCheck;
2281 template<
typename SrcTreeT,
typename LeafManagerT,
typename AuxDataT>
2283 const LeafManagerT& leafs,
double iso,
bool extraCheck)
2284 : mLeafManager(leafs)
2286 , mSourceAccessor(mSourceTree)
2288 , mEdgeAccessor(*mEdgeTree)
2289 , mAuxTree(new
AuxTreeT(AuxDataT(0)))
2290 , mAuxAccessor(*mAuxTree)
2292 , mExtraCheck(extraCheck)
2296 template<
typename SrcTreeT,
typename LeafManagerT,
typename AuxDataT>
2298 : mLeafManager(rhs.mLeafManager)
2299 , mSourceTree(rhs.mSourceTree)
2300 , mSourceAccessor(mSourceTree)
2302 , mEdgeAccessor(*mEdgeTree)
2303 , mAuxTree(new
AuxTreeT(AuxDataT(0)))
2304 , mAuxAccessor(*mAuxTree)
2305 , mIsovalue(rhs.mIsovalue)
2306 , mExtraCheck(rhs.mExtraCheck)
2311 template<
typename SrcTreeT,
typename LeafManagerT,
typename AuxDataT>
2316 tbb::parallel_reduce(mLeafManager.getRange(), *
this);
2318 (*this)(mLeafManager.getRange());
2323 template<
typename SrcTreeT,
typename LeafManagerT,
typename AuxDataT>
2327 typename LeafManagerT::TreeType::LeafNodeType::ValueOnCIter iter;
2334 for (
size_t n = range.begin(); n != range.end(); ++n) {
2335 for (iter = mLeafManager.leaf(n).cbeginValueOn(); iter; ++iter) {
2337 ijk = iter.getCoord();
2339 thisInside = mSourceAccessor.getValue(ijk) < mIsovalue;
2340 edgeFlags = edgeCheck(ijk, thisInside);
2342 if (edgeFlags != 0) {
2343 edgeFlags |= int(thisInside);
2344 mEdgeAccessor.setValue(ijk,
char(edgeFlags));
2349 for (
size_t n = range.begin(); n != range.end(); ++n) {
2350 for (iter = mLeafManager.leaf(n).cbeginValueOn(); iter; ++iter) {
2352 ijk = iter.getCoord();
2354 thisInside = mSourceAccessor.getValue(ijk) < mIsovalue;
2355 edgeFlags = edgeCheck(ijk, thisInside);
2357 if (edgeFlags != 0) {
2358 edgeFlags |= int(thisInside);
2359 mEdgeAccessor.setValue(ijk,
char(edgeFlags));
2363 if (!mSourceAccessor.probeValue(ijk, val)) {
2364 thisInside = val < mIsovalue;
2365 edgeFlags = edgeCheck(ijk, thisInside);
2367 if (edgeFlags != 0) {
2368 edgeFlags |= int(thisInside);
2369 mEdgeAccessor.setValue(ijk,
char(edgeFlags));
2375 if (!mSourceAccessor.probeValue(ijk, val)) {
2376 thisInside = val < mIsovalue;
2377 edgeFlags = edgeCheck(ijk, thisInside);
2379 if (edgeFlags != 0) {
2380 edgeFlags |= int(thisInside);
2381 mEdgeAccessor.setValue(ijk,
char(edgeFlags));
2387 if (!mSourceAccessor.probeValue(ijk, val)) {
2388 thisInside = val < mIsovalue;
2389 edgeFlags = edgeCheck(ijk, thisInside);
2391 if (edgeFlags != 0) {
2392 edgeFlags |= int(thisInside);
2393 mEdgeAccessor.setValue(ijk,
char(edgeFlags));
2401 template<
typename SrcTreeT,
typename LeafManagerT,
typename AuxDataT>
2409 n_ijk = ijk; ++n_ijk[0];
2410 bool otherInside = (mSourceAccessor.getValue(n_ijk) < mIsovalue);
2411 if (otherInside != thisInside) {
2415 mAuxAccessor.setActiveState(ijk,
true);
2417 coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2];
2418 mAuxAccessor.setActiveState(coord,
true);
2420 coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2]-1;
2421 mAuxAccessor.setActiveState(coord,
true);
2423 coord[0] = ijk[0]; coord[1] = ijk[1]; coord[2] = ijk[2]-1;
2424 mAuxAccessor.setActiveState(coord,
true);
2428 n_ijk[0] = ijk[0]; ++n_ijk[1];
2429 otherInside = (mSourceAccessor.getValue(n_ijk) < mIsovalue);
2430 if (otherInside != thisInside) {
2434 mAuxAccessor.setActiveState(ijk,
true);
2436 coord[0] = ijk[0]; coord[1] = ijk[1]; coord[2] = ijk[2]-1;
2437 mAuxAccessor.setActiveState(coord,
true);
2439 coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2];
2440 mAuxAccessor.setActiveState(coord,
true);
2442 coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2]-1;
2443 mAuxAccessor.setActiveState(coord,
true);
2447 n_ijk[1] = ijk[1]; ++n_ijk[2];
2448 otherInside = (mSourceAccessor.getValue(n_ijk) < mIsovalue);
2449 if (otherInside != thisInside) {
2453 mAuxAccessor.setActiveState(ijk,
true);
2455 coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2];
2456 mAuxAccessor.setActiveState(coord,
true);
2458 coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2];
2459 mAuxAccessor.setActiveState(coord,
true);
2461 coord[0] = ijk[0]-1; coord[1] = ijk[1]-1; coord[2] = ijk[2];
2462 mAuxAccessor.setActiveState(coord,
true);
2473 template <
class DistTreeT>
2477 typedef typename DistTreeT::template ValueConverter<bool>::Type
BoolTreeT;
2478 typedef typename DistTreeT::template ValueConverter<int>::Type
IntTreeT;
2490 void operator()(
const tbb::blocked_range<size_t>&)
const;
2501 template <
class DistTreeT>
2504 : mSeamMaskLeafs(seamMaskLeafs)
2505 , mTopologyMaskTree(topologyMaskTree)
2506 , mTopologyMaskAcc(mTopologyMaskTree)
2512 template <
class DistTreeT>
2514 : mSeamMaskLeafs(rhs.mSeamMaskLeafs)
2515 , mTopologyMaskTree(rhs.mTopologyMaskTree)
2516 , mTopologyMaskAcc(mTopologyMaskTree)
2517 , mAuxTree(rhs.mAuxTree)
2522 template <
class DistTreeT>
2526 tbb::parallel_for(mSeamMaskLeafs.getRange(), *
this);
2529 template <
class DistTreeT>
2533 (*this)(mSeamMaskLeafs.getRange());
2536 template <
class DistTreeT>
2540 typedef typename BoolTreeT::LeafNodeType::ValueOnIter ValueOnIterT;
2542 for (
size_t leafIdx = range.begin(); leafIdx != range.end(); ++leafIdx) {
2543 ValueOnIterT it = mSeamMaskLeafs[leafIdx]->beginValueOn();
2545 ijk = it.getCoord();
2546 if (!mTopologyMaskAcc.isValueOn(ijk)) {
2549 bool turnOff =
true;
2550 for (
size_t n = 0; n < 6; ++n) {
2552 if (!mAuxTree.isValueOn(n_ijk) && mTopologyMaskAcc.isValueOn(n_ijk)) {
2557 if (turnOff) it.setValueOff();
2566 template <
class DistTreeT>
2570 typedef typename DistTreeT::template ValueConverter<bool>::Type
BoolTreeT;
2571 typedef typename DistTreeT::template ValueConverter<int>::Type
IntTreeT;
2585 void runSerial(
const size_t iterations);
2587 void operator()(
const tbb::blocked_range<size_t>&)
const;
2596 bool pointInAABB(
const Vec3s& p,
const Vec3s& bmin,
const Vec3s& bmax)
const
2598 for (
int i = 0; i < 3; ++i) {
2599 if (p[i] < bmin[i] || p[i] > bmax[i]) {
2609 template <
class DistTreeT>
2617 , mEdgeMaskTree(edgeMaskTree)
2624 template <
class DistTreeT>
2626 : mLeafs(rhs.mLeafs)
2627 , mEdgeMaskTree(rhs.mEdgeMaskTree)
2628 , mAuxTree(rhs.mAuxTree)
2629 , mPoints(rhs.mPoints)
2630 , mTransform(rhs.mTransform)
2634 template <
class DistTreeT>
2638 for (
size_t i = 0; i < iterations; ++i) {
2639 tbb::parallel_for(mLeafs.getRange(), *
this);
2643 template <
class DistTreeT>
2647 for (
size_t i = 0; i < iterations; ++i) {
2648 (*this)(mLeafs.getRange());
2652 template <
class DistTreeT>
2656 typedef typename BoolTreeT::LeafNodeType::ValueOnIter ValueOnIterT;
2660 float dx = float(mTransform.voxelSize()[0]);
2665 for (
size_t leafIdx = range.begin(); leafIdx != range.end(); ++leafIdx) {
2666 typename BoolTreeT::LeafNodeType::ValueOnIter valueIt = mLeafs[leafIdx]->beginValueOn();
2667 for ( ; valueIt; ++valueIt) {
2669 ijk = valueIt.getCoord();
2674 for (
int n = 0; n < 26; ++n) {
2677 avg += mPoints[auxAcc.
getValue(nijk)];
2683 avg *= (1.0 / float(count));
2689 bool inCell =
false;
2690 for (
int i = 0; i < 10; ++i) {
2692 inCell = pointInAABB(avg, bmin, bmax);
2702 if (inCell) ptn = avg;
2731 template<
class CharAccessorT,
typename AuxAccessorT>
2736 : mEdgeAcc(edgeAcc), mAuxAcc(auxAcc) {}
2741 mEdgeAcc.setValue(ijk, edgeFlags |
XEDGE);
2743 mAuxAcc.setActiveState(ijk,
true);
2746 coord[1] = ijk[1]-1;
2747 mAuxAcc.setActiveState(coord,
true);
2749 coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2]-1;
2750 mAuxAcc.setActiveState(coord,
true);
2752 coord[0] = ijk[0]; coord[1] = ijk[1]; coord[2] = ijk[2]-1;
2753 mAuxAcc.setActiveState(coord,
true);
2758 mEdgeAcc.setValue(ijk, edgeFlags |
YEDGE);
2760 mAuxAcc.setActiveState(ijk,
true);
2763 coord[2] = ijk[2]-1;
2764 mAuxAcc.setActiveState(coord,
true);
2766 coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2];
2767 mAuxAcc.setActiveState(coord,
true);
2769 coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2]-1;
2770 mAuxAcc.setActiveState(coord,
true);
2775 mEdgeAcc.setValue(ijk, edgeFlags |
ZEDGE);
2777 mAuxAcc.setActiveState(ijk,
true);
2780 coord[1] = ijk[1]-1;
2781 mAuxAcc.setActiveState(coord,
true);
2783 coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2];
2784 mAuxAcc.setActiveState(coord,
true);
2786 coord[0] = ijk[0]-1; coord[1] = ijk[1]-1; coord[2] = ijk[2];
2787 mAuxAcc.setActiveState(coord,
true);
2791 CharAccessorT& mEdgeAcc;
2792 AuxAccessorT& mAuxAcc;
2799 template<
class DistTreeT,
class AuxTreeT,
class CharTreeT>
2802 const DistTreeT& distTree, CharTreeT& edgeTree, AuxTreeT& auxTree,
2809 typename DistTreeT::ValueType isoValue =
typename DistTreeT::ValueType(iso);
2811 DistAccessorT distAcc(distTree);
2812 CharTreeAccessorT edgeAcc(edgeTree);
2813 AuxTreeAccessorT auxAcc(auxTree);
2818 typename DistTreeT::ValueType value;
2820 bool processTileFace;
2822 typename DistTreeT::ValueOnCIter tileIter(distTree);
2823 tileIter.setMaxDepth(DistTreeT::ValueOnCIter::LEAF_DEPTH - 1);
2825 for ( ; tileIter; ++tileIter) {
2826 tileIter.getBoundingBox(bbox);
2828 const bool thisInside = (distAcc.getValue(bbox.
min()) < isoValue);
2829 const int thisDepth = distAcc.getValueDepth(bbox.
min());
2838 processTileFace =
true;
2839 if (thisDepth >= distAcc.getValueDepth(nijk)) {
2840 processTileFace = thisInside != (distAcc.getValue(nijk) < isoValue);
2843 if (processTileFace) {
2844 for (ijk[1] = bbox.
min()[1]; ijk[1] <= bbox.
max()[1]; ++ijk[1]) {
2846 for (ijk[2] = bbox.
min()[2]; ijk[2] <= bbox.
max()[2]; ++ijk[2]) {
2848 if ((distAcc.getValue(nijk) < isoValue) != thisInside) {
2849 auxData.
setXEdge(edgeAcc.getValue(ijk) | char(thisInside), ijk);
2858 processTileFace =
true;
2859 if (thisDepth >= distAcc.getValueDepth(ijk)) {
2860 processTileFace = !distAcc.probeValue(ijk, value) && thisInside != (value < isoValue);
2863 if (processTileFace) {
2864 for (ijk[1] = bbox.
min()[1]; ijk[1] <= bbox.
max()[1]; ++ijk[1]) {
2865 for (ijk[2] = bbox.
min()[2]; ijk[2] <= bbox.
max()[2]; ++ijk[2]) {
2866 if (!distAcc.probeValue(ijk, value) && (value < isoValue) != thisInside) {
2867 auxData.
setXEdge(edgeAcc.getValue(ijk) | char(!thisInside), ijk);
2880 processTileFace =
true;
2881 if (thisDepth >= distAcc.getValueDepth(nijk)) {
2882 processTileFace = thisInside != (distAcc.getValue(nijk) < isoValue);
2885 if (processTileFace) {
2886 for (ijk[0] = bbox.
min()[0]; ijk[0] <= bbox.
max()[0]; ++ijk[0]) {
2888 for (ijk[2] = bbox.
min()[2]; ijk[2] <= bbox.
max()[2]; ++ijk[2]) {
2891 if ((distAcc.getValue(nijk) < isoValue) != thisInside) {
2892 auxData.
setYEdge(edgeAcc.getValue(ijk) | char(thisInside), ijk);
2902 processTileFace =
true;
2903 if (thisDepth >= distAcc.getValueDepth(ijk)) {
2904 processTileFace = !distAcc.probeValue(ijk, value) && thisInside != (value < isoValue);
2907 if (processTileFace) {
2908 for (ijk[0] = bbox.
min()[0]; ijk[0] <= bbox.
max()[0]; ++ijk[0]) {
2909 for (ijk[2] = bbox.
min()[2]; ijk[2] <= bbox.
max()[2]; ++ijk[2]) {
2911 if (!distAcc.probeValue(ijk, value) && (value < isoValue) != thisInside) {
2912 auxData.
setYEdge(edgeAcc.getValue(ijk) | char(!thisInside), ijk);
2925 processTileFace =
true;
2926 if (thisDepth >= distAcc.getValueDepth(nijk)) {
2927 processTileFace = thisInside != (distAcc.getValue(nijk) < isoValue);
2930 if (processTileFace) {
2931 for (ijk[0] = bbox.
min()[0]; ijk[0] <= bbox.
max()[0]; ++ijk[0]) {
2933 for (ijk[1] = bbox.
min()[1]; ijk[1] <= bbox.
max()[1]; ++ijk[1]) {
2936 if ((distAcc.getValue(nijk) < isoValue) != thisInside) {
2937 auxData.
setZEdge(edgeAcc.getValue(ijk) | char(thisInside), ijk);
2946 processTileFace =
true;
2947 if (thisDepth >= distAcc.getValueDepth(ijk)) {
2948 processTileFace = !distAcc.probeValue(ijk, value) && thisInside != (value < isoValue);
2951 if (processTileFace) {
2952 for (ijk[0] = bbox.
min()[0]; ijk[0] <= bbox.
max()[0]; ++ijk[0]) {
2953 for (ijk[1] = bbox.
min()[1]; ijk[1] <= bbox.
max()[1]; ++ijk[1]) {
2955 if (!distAcc.probeValue(ijk, value) && (value < isoValue) != thisInside) {
2956 auxData.
setZEdge(edgeAcc.getValue(ijk) | char(!thisInside), ijk);
2973 : mPointsIn(pointsIn) , mPointsOut(pointsOut)
2979 for (
size_t n = range.begin(); n < range.end(); ++n) {
2980 mPointsOut[n] = mPointsIn[n];
2986 std::vector<Vec3s>& mPointsOut;
2998 , mPolygonPoolListSize(0)
2999 , mIsovalue(isovalue)
3000 , mPrimAdaptivity(adaptivity)
3001 , mSecAdaptivity(0.0)
3003 , mSurfaceMaskGrid(
GridBase::ConstPtr())
3004 , mAdaptivityGrid(
GridBase::ConstPtr())
3005 , mAdaptivityMaskTree(
TreeBase::ConstPtr())
3007 , mRefTopologyMaskTree(
TreeBase::Ptr())
3009 , mSmoothSeams(false)
3010 , mInvertSurfaceMask(false)
3022 inline const size_t&
3025 return mPointListSize;
3043 inline const size_t&
3046 return mPolygonPoolListSize;
3054 mSecAdaptivity = secAdaptivity;
3059 mSmoothSeams = smoothSeams;
3065 mSurfaceMaskGrid = mask;
3066 mInvertSurfaceMask = invertMask;
3072 mAdaptivityGrid = grid;
3078 mAdaptivityMaskTree = tree;
3084 mPartitions =
std::max(partitions,
unsigned(1));
3085 mActivePart =
std::min(activePart, mPartitions-1);
3091 template<
typename Gr
idT>
3095 typedef typename GridT::TreeType DistTreeT;
3096 typedef typename DistTreeT::ValueType DistValueT;
3098 typedef typename DistTreeT::template ValueConverter<char>::Type CharTreeT;
3099 typedef typename DistTreeT::template ValueConverter<int>::Type IntTreeT;
3100 typedef typename DistTreeT::template ValueConverter<Vec3s>::Type Vec3sTreeT;
3101 typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
3104 typedef typename DistTreeT::template ValueConverter<float>::Type FloatTreeT;
3107 const openvdb::math::Transform& transform = distGrid.
transform();
3108 const DistTreeT& distTree = distGrid.tree();
3109 typename CharTreeT::Ptr edgeTree;
3110 typename IntTreeT::Ptr auxTree;
3111 typename BoolTreeT::Ptr pointMask;
3113 const bool nonLevelSetGrid = distGrid.getGridClass() !=
GRID_LEVEL_SET;
3114 const bool noAdaptivity = mPrimAdaptivity < 1e-6 && mSecAdaptivity < 1e-6;
3116 const bool extraSignCheck = nonLevelSetGrid ||
3117 (std::abs(mIsovalue -
double(distGrid.background())) < (1.5 * transform.voxelSize()[0]));
3119 const BoolGridT * surfaceMask = NULL;
3120 if (mSurfaceMaskGrid) {
3121 surfaceMask =
static_cast<const BoolGridT*
>(mSurfaceMaskGrid.get());
3124 const FloatGridT * adaptivityField = NULL;
3125 if (mAdaptivityGrid) {
3126 adaptivityField =
static_cast<const FloatGridT*
>(mAdaptivityGrid.get());
3131 if (surfaceMask || mPartitions > 1) {
3133 BoolTreeT valueMask(
false);
3141 *surfaceMask, leafs, transform, mInvertSurfaceMask);
3143 valueMask.merge(masking.
tree());
3147 if (mPartitions > 1) {
3150 valueMask.pruneInactive();
3159 valueMask.merge(partitioner.
tree());
3165 op(distTree, leafs, mIsovalue, extraSignCheck);
3173 if (!noAdaptivity) {
3174 pointMask =
typename BoolTreeT::Ptr(
new BoolTreeT(
false));
3175 pointMask->topologyUnion(*auxTree);
3188 op(distTree, leafs, mIsovalue, extraSignCheck);
3208 if (nonLevelSetGrid) {
3217 const GridT* refGrid =
static_cast<const GridT*
>(mRefGrid.get());
3218 if (refGrid && refGrid->activeVoxelCount() > 0) {
3224 if (!mRefEdgeTree && !mRefTopologyMaskTree) {
3227 *refData.
mDistTree, leafs, mIsovalue, extraSignCheck);
3230 mRefTopologyMaskTree = op.
auxTree();
3231 mSeamPointTree =
typename Vec3sTreeT::Ptr(
new Vec3sTreeT(
Vec3s(0.0)));
3234 if (mRefEdgeTree && mRefTopologyMaskTree) {
3235 refData.
mEdgeTree =
static_cast<CharTreeT*
>(mRefEdgeTree.get());
3236 refData.
mTopologyMaskTree =
static_cast<BoolTreeT*
>(mRefTopologyMaskTree.get());
3237 refData.
mSeamPointTree =
static_cast<Vec3sTreeT*
>(mSeamPointTree.get());
3238 refData.
mSeamMaskTree =
typename BoolTreeT::Ptr(
new BoolTreeT(
false));
3244 BoolTreeT edgeMaskTree(0.0);
3263 if (mAdaptivityMaskTree) {
3265 const BoolTreeT* refGrid =
static_cast<const BoolTreeT*
>(
3266 mAdaptivityMaskTree.get());
3275 std::vector<size_t> regions(auxLeafs.
size(), 0);
3283 regions, DistValueT(mIsovalue), DistValueT(mPrimAdaptivity), pointMask.get());
3287 if (adaptivityField) {
3288 merge.setSpatialAdaptivity(transform, *adaptivityField);
3296 for (
size_t n = 0, N = regions.size(); n < N; ++n) {
3298 regions[n] = mPointListSize;
3299 mPointListSize += tmp;
3305 for (
size_t n = 0, N = auxLeafs.
size(); n < N; ++n) {
3306 acc.
touchLeaf(auxLeafs[n]->getOrigin());
3314 pointGen(distTree, auxLeafs, regions, transform, mPoints, mIsovalue);
3320 if (mSmoothSeams && refData.
isValid()) {
3332 typename BoolTreeT::LeafIter leafIt = refData.
mSeamMaskTree->beginLeaf();
3333 for ( ; leafIt; ++leafIt) {
3334 typename BoolTreeT::LeafNodeType::ValueOnIter valueIt = leafIt->beginValueOn();
3335 for ( ; valueIt; ++valueIt) {
3336 ijk = valueIt.getCoord();
3337 const int idx = auxAcc.
getValue(ijk);
3338 if (idx != 0 && !ptnAcc.
isValueOn(ijk)) {
3339 ptnAcc.
setValue(ijk, mPoints[idx]);
3348 mPolygonPoolListSize = edgeLeafs.
size();
3349 mPolygons.reset(
new PolygonPool[mPolygonPoolListSize]);
3353 meshGen(edgeLeafs, *auxTree, mPolygons);
3358 meshGen(edgeLeafs, *auxTree, mPolygons);
3369 template<
typename Gr
idType>
3372 const GridType& grid,
3373 std::vector<Vec3s>& points,
3374 std::vector<Vec3I>& triangles,
3375 std::vector<Vec4I>& quads,
3388 tbb::parallel_for(tbb::blocked_range<size_t>(0, points.size()), ptnCpy);
3395 size_t numQuads = 0, numTriangles = 0;
3397 openvdb::tools::PolygonPool& polygons = polygonPoolList[n];
3398 numTriangles += polygons.numTriangles();
3399 numQuads += polygons.numQuads();
3403 triangles.resize(numTriangles);
3405 quads.resize(numQuads);
3409 size_t qIdx = 0, tIdx = 0;
3411 openvdb::tools::PolygonPool& polygons = polygonPoolList[n];
3413 for (
size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
3414 quads[qIdx++] = polygons.quad(i);
3417 for (
size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
3418 triangles[tIdx++] = polygons.triangle(i);
3424 template<
typename Gr
idType>
3427 const GridType& grid,
3428 std::vector<Vec3s>& points,
3429 std::vector<Vec4I>& quads,
3432 std::vector<Vec3I> triangles(0);
3433 volumeToMesh(grid,points, triangles, quads, isovalue, 0.0);
3440 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED