31 #ifndef OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
32 #define OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
34 #include <openvdb/Types.h>
35 #include <openvdb/math/FiniteDifference.h>
36 #include <openvdb/math/Operators.h>
37 #include <openvdb/math/Proximity.h>
38 #include <openvdb/tools/LevelSetUtil.h>
39 #include <openvdb/tools/Morphology.h>
40 #include <openvdb/util/NullInterrupter.h>
41 #include <openvdb/util/Util.h>
43 #include <tbb/blocked_range.h>
44 #include <tbb/parallel_for.h>
45 #include <tbb/parallel_reduce.h>
61 template<
typename DistGr
idT,
typename InterruptT = util::NullInterrupter>
69 typedef typename DistTreeT::template ValueConverter<Int32>::Type
IndexTreeT;
71 typedef typename DistTreeT::template ValueConverter<bool>::Type
StencilTreeT;
75 MeshToVolume(openvdb::math::Transform::Ptr&,
int conversionFlags = 0,
76 InterruptT *interrupter = NULL,
int signSweeps = 1);
89 void convertToLevelSet(
const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList,
101 void convertToUnsignedDistanceField(
const std::vector<Vec3s>& pointList,
102 const std::vector<Vec4I>& polygonList,
DistValueT exBandWidth);
126 void convert(
const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList,
127 DistValueT exBandWidth, DistValueT inBandWidth = DistValueT(0.0));
135 void doConvert(
const std::vector<Vec3s>&,
const std::vector<Vec4I>&,
136 DistValueT exBandWidth, DistValueT inBandWidth,
bool unsignedDistField =
false);
138 openvdb::math::Transform::Ptr mTransform;
139 int mConversionFlags, mSignSweeps;
141 typename DistGridT::Ptr mDistGrid;
142 typename IndexGridT::Ptr mIndexGrid;
143 typename StencilGridT::Ptr mIntersectingVoxelsGrid;
145 InterruptT *mInterrupter;
157 template<
typename DistTreeT,
typename IndexTreeT>
159 combine(DistTreeT& lhsDist, IndexTreeT& lhsIndex, DistTreeT& rhsDist, IndexTreeT& rhsIndex)
161 typedef typename DistTreeT::ValueType DistValueT;
165 typename DistTreeT::LeafCIter iter = rhsDist.cbeginLeaf();
170 for ( ; iter; ++iter) {
171 typename DistTreeT::LeafNodeType::ValueOnCIter it = iter->cbeginValueOn();
176 rhsValue = it.getValue();
177 DistValueT& lhsValue =
const_cast<DistValueT&
>(lhsDistAccessor.
getValue(ijk));
179 if (-rhsValue < std::abs(lhsValue)) {
198 template<
typename DistTreeT,
typename InterruptT = util::NullInterrupter>
206 typedef typename DistTreeT::template ValueConverter<Int32>::Type
IndexTreeT;
208 typedef typename DistTreeT::template ValueConverter<bool>::Type
StencilTreeT;
213 const std::vector<Vec4I>& polygonList, InterruptT *interrupter = NULL);
221 void operator() (
const tbb::blocked_range<size_t> &range);
233 bool evalVoxel(
const Coord& ijk,
const Int32 polyIdx);
235 std::vector<Vec3s>
const *
const mPointList;
236 std::vector<Vec4I>
const *
const mPolygonList;
238 DistTreeT mSqrDistTree;
239 DistAccessorT mSqrDistAccessor;
241 IndexTreeT mPrimIndexTree;
242 IndexAccessorT mPrimIndexAccessor;
244 StencilTreeT mIntersectionTree;
245 StencilAccessorT mIntersectionAccessor;
248 IndexTreeT mLastPrimTree;
249 IndexAccessorT mLastPrimAccessor;
251 InterruptT *mInterrupter;
255 template<
typename DistTreeT,
typename InterruptT>
259 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPolygonList->size()), *
this);
262 template<
typename DistTreeT,
typename InterruptT>
266 (*this)(tbb::blocked_range<size_t>(0, mPolygonList->size()));
269 template<
typename DistTreeT,
typename InterruptT>
271 const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList,
272 InterruptT *interrupter)
273 : mPointList(&pointList)
274 , mPolygonList(&polygonList)
276 , mSqrDistAccessor(mSqrDistTree)
278 , mPrimIndexAccessor(mPrimIndexTree)
279 , mIntersectionTree(false)
280 , mIntersectionAccessor(mIntersectionTree)
282 , mLastPrimAccessor(mLastPrimTree)
283 , mInterrupter(interrupter)
287 template<
typename DistTreeT,
typename InterruptT>
290 : mPointList(rhs.mPointList)
291 , mPolygonList(rhs.mPolygonList)
293 , mSqrDistAccessor(mSqrDistTree)
295 , mPrimIndexAccessor(mPrimIndexTree)
296 , mIntersectionTree(false)
297 , mIntersectionAccessor(mIntersectionTree)
299 , mLastPrimAccessor(mLastPrimTree)
300 , mInterrupter(rhs.mInterrupter)
304 template<
typename DistTreeT,
typename InterruptT>
309 double edge_max = std::abs(v1[0] - v0[0]);
310 edge_max =
std::max(edge_max, std::abs(v1[1] - v0[1]));
311 edge_max =
std::max(edge_max, std::abs(v1[2] - v0[2]));
312 edge_max =
std::max(edge_max, std::abs(v0[0] - v2[0]));
313 edge_max =
std::max(edge_max, std::abs(v0[1] - v2[1]));
314 edge_max =
std::max(edge_max, std::abs(v0[2] - v2[2]));
315 return edge_max < 200.0;
318 template<
typename DistTreeT,
typename InterruptT>
322 std::deque<Coord> coordList;
327 for (
size_t n = range.begin(); n < range.end(); ++n) {
329 if (mInterrupter && mInterrupter->wasInterrupted()) {
330 tbb::task::self().cancel_group_execution();
334 const Int32 primIdx = n;
335 const Vec4I verts = (*mPolygonList)[n];
337 Vec3d p0((*mPointList)[verts[0]]);
338 Vec3d p1((*mPointList)[verts[1]]);
339 Vec3d p2((*mPointList)[verts[2]]);
341 if (shortEdge(p0, p1, p2)) {
346 evalVoxel(ijk, primIdx);
347 coordList.push_back(ijk);
351 evalVoxel(ijk, primIdx);
352 coordList.push_back(ijk);
356 evalVoxel(ijk, primIdx);
357 coordList.push_back(ijk);
360 Vec3d p3((*mPointList)[verts[3]]);
362 evalVoxel(ijk, primIdx);
363 coordList.push_back(ijk);
366 while (!coordList.empty()) {
367 if (mInterrupter && mInterrupter->wasInterrupted()) {
371 ijk = coordList.back();
372 coordList.pop_back();
374 mIntersectionAccessor.setActiveState(ijk,
true);
376 for (
Int32 i = 0; i < 26; ++i) {
379 if (primIdx != mLastPrimAccessor.getValue(n_ijk)) {
380 mLastPrimAccessor.setValue(n_ijk, n);
381 if(evalVoxel(n_ijk, n)) coordList.push_back(n_ijk);
389 evalVoxel(ijk, primIdx);
391 mLastPrimAccessor.setValue(ijk, primIdx);
394 while (!auxTree.empty()) {
396 if (mInterrupter && mInterrupter->wasInterrupted()) {
400 typename StencilTreeT::LeafIter leafIter = auxTree.beginLeaf();
401 for (; leafIter; leafIter.next()) {
403 if (mInterrupter && mInterrupter->wasInterrupted()) {
407 typename StencilTreeT::LeafNodeType::ValueOnIter iter = leafIter->beginValueOn();
408 for (; iter; iter.next()) {
409 ijk = iter.getCoord();
412 mIntersectionAccessor.setActiveState(ijk,
true);
414 for (
Int32 i = 0; i < 26; ++i) {
417 if (primIdx != mLastPrimAccessor.getValue(n_ijk)) {
418 mLastPrimAccessor.setValue(n_ijk, n);
425 auxTree.pruneInactive();
431 template<
typename DistTreeT,
typename InterruptT>
435 Vec3d voxelCenter(ijk[0], ijk[1], ijk[2]);
436 Vec4I verts = (*mPolygonList)[polyIdx];
439 Vec3d p0((*mPointList)[verts[0]]);
440 Vec3d p1((*mPointList)[verts[1]]);
441 Vec3d p2((*mPointList)[verts[2]]);
447 p1 =
Vec3d((*mPointList)[verts[3]]);
450 if (secondDist < dist) dist = secondDist;
453 const DistValueT tmp(dist);
454 if (tmp < std::abs(mSqrDistAccessor.getValue(ijk))) {
455 mSqrDistAccessor.setValue(ijk, -tmp);
456 mPrimIndexAccessor.setValue(ijk, polyIdx);
459 return (dist < 0.86602540378443861);
462 template<
typename DistTreeT,
typename InterruptT>
466 typename DistTreeT::LeafCIter iter = rhs.mSqrDistTree.cbeginLeaf();
470 for ( ; iter; ++iter) {
471 typename DistTreeT::LeafNodeType::ValueOnCIter it = iter->cbeginValueOn();
476 rhsDist = it.getValue();
477 DistValueT lhsDist = mSqrDistAccessor.getValue(ijk);
479 if (-rhsDist < std::abs(lhsDist)) {
480 mSqrDistAccessor.setValue(ijk, rhsDist);
481 mPrimIndexAccessor.setValue(ijk, rhs.mPrimIndexAccessor.getValue(ijk));
486 mIntersectionTree.merge(rhs.mIntersectionTree);
496 template<
typename DistTreeT,
typename InterruptT = util::NullInterrupter>
504 typedef typename DistTreeT::template ValueConverter<bool>::Type
StencilTreeT;
515 void operator()(
const tbb::blocked_range<int> &range)
const;
520 int sparseScan(
int slice)
const;
522 DistTreeT& mDistTree;
531 std::vector<Index> mStepSize;
533 InterruptT *mInterrupter;
536 template<
typename DistTreeT,
typename InterruptT>
540 tbb::parallel_for(tbb::blocked_range<int>(mBBox.min()[0], mBBox.max()[0]+1), *
this);
543 template<
typename DistTreeT,
typename InterruptT>
547 (*this)(tbb::blocked_range<int>(mBBox.min()[0], mBBox.max()[0]+1));
550 template<
typename DistTreeT,
typename InterruptT>
552 DistTreeT& distTree,
const StencilTreeT& intersectionTree, InterruptT *interrupter)
553 : mDistTree(distTree)
554 , mDistAccessor(mDistTree)
555 , mIntersectionTree(intersectionTree)
556 , mIntersectionAccessor(mIntersectionTree)
559 , mInterrupter(interrupter)
562 std::vector<Index> dims;
563 mDistTree.getNodeLog2Dims(dims);
565 mStepSize.resize(dims.size()+1, 1);
567 for (
int idx = static_cast<int>(dims.size()) - 1; idx > -1; --idx) {
568 exponent += dims[idx];
569 mStepSize[idx] = 1 << exponent;
572 mDistTree.evalLeafBoundingBox(mBBox);
575 const int tileDim = mStepSize[0];
577 for (
size_t i = 0; i < 3; ++i) {
580 double diff = std::abs(
double(mBBox.
min()[i])) / double(tileDim);
582 if (mBBox.
min()[i] <= tileDim) {
583 n = int(std::ceil(diff));
584 mBBox.
min()[i] = - n * tileDim;
586 n = int(std::floor(diff));
587 mBBox.
min()[i] = n * tileDim;
590 n = int(std::ceil(std::abs(
double(mBBox.
max()[i] - mBBox.
min()[i])) / double(tileDim)));
591 mBBox.
max()[i] = mBBox.
min()[i] + n * tileDim;
595 template<
typename DistTreeT,
typename InterruptT>
598 : mDistTree(rhs.mDistTree)
599 , mDistAccessor(mDistTree)
600 , mIntersectionTree(rhs.mIntersectionTree)
601 , mIntersectionAccessor(mIntersectionTree)
603 , mStepSize(rhs.mStepSize)
604 , mInterrupter(rhs.mInterrupter)
608 template<
typename DistTreeT,
typename InterruptT>
614 for (
int n = range.begin(); n < range.end(); n += iStep) {
616 if (mInterrupter && mInterrupter->wasInterrupted()) {
617 tbb::task::self().cancel_group_execution();
621 iStep = sparseScan(n);
625 template<
typename DistTreeT,
typename InterruptT>
629 bool lastVoxelWasOut =
true;
632 Coord ijk(slice, mBBox.min()[1], mBBox.min()[2]);
633 Coord step(mStepSize[mDistAccessor.getValueDepth(ijk) + 1]);
636 for (ijk[1] = mBBox.
min()[1]; ijk[1] <= mBBox.
max()[1]; ijk[1] += step[1]) {
638 if (mInterrupter && mInterrupter->wasInterrupted()) {
642 step[1] = mStepSize[mDistAccessor.getValueDepth(ijk) + 1];
643 step[0] =
std::min(step[0], step[1]);
645 for (ijk[2] = mBBox.
min()[2]; ijk[2] <= mBBox.
max()[2]; ijk[2] += step[2]) {
647 step[2] = mStepSize[mDistAccessor.getValueDepth(ijk) + 1];
648 step[1] =
std::min(step[1], step[2]);
649 step[0] =
std::min(step[0], step[2]);
652 if (mDistAccessor.isValueOn(ijk)) {
655 if (mIntersectionAccessor.isValueOn(ijk)) {
657 lastVoxelWasOut =
false;
660 }
else if (lastVoxelWasOut) {
662 DistValueT& val =
const_cast<DistValueT&
>(mDistAccessor.getValue(ijk));
668 for (
Int32 n = 3; n < 6; n += 2) {
671 if (mDistAccessor.probeValue(n_ijk, val) && val > 0) {
672 lastVoxelWasOut =
true;
677 if (lastVoxelWasOut) {
679 DistValueT& v =
const_cast<DistValueT&
>(mDistAccessor.getValue(ijk));
682 const int tmp_k = ijk[2];
685 for (--ijk[2]; ijk[2] >= last_k; --ijk[2]) {
686 if (mIntersectionAccessor.isValueOn(ijk))
break;
687 DistValueT& v =
const_cast<DistValueT&
>(mDistAccessor.getValue(ijk));
688 if(v < DistValueT(0.0)) v = -v;
714 template<
typename DistTreeT>
722 typedef typename DistTreeT::template ValueConverter<Int32>::Type
IndexTreeT;
724 typedef typename DistTreeT::template ValueConverter<bool>::Type
StencilTreeT;
730 const std::vector<Vec3s>& pointList,
731 const std::vector<Vec4I>& polygonList,
743 void operator()(
const tbb::blocked_range<size_t>&)
const;
748 void evalVoxel(
const Coord& ijk)
const;
749 Vec3d getClosestPointDir(
const Coord& ijk)
const;
751 std::vector<Vec3s>
const *
const mPointList;
752 std::vector<Vec4I>
const *
const mPolygonList;
754 DistTreeT& mDistTree;
765 template<
typename DistTreeT>
769 tbb::parallel_for(
mLeafs.getRange(), *
this);
772 template<
typename DistTreeT>
776 (*this)(
mLeafs.getRange());
779 template<
typename DistTreeT>
781 const std::vector<Vec3s>& pointList,
782 const std::vector<Vec4I>& polygonList,
787 : mPointList(&pointList)
788 , mPolygonList(&polygonList)
789 , mDistTree(distTree)
790 , mDistAccessor(mDistTree)
791 , mIndexTree(indexTree)
792 , mIndexAccessor(mIndexTree)
793 , mIntersectionTree(intersectionTree)
794 , mIntersectionAccessor(mIntersectionTree)
799 template<
typename DistTreeT>
802 : mPointList(rhs.mPointList)
803 , mPolygonList(rhs.mPolygonList)
804 , mDistTree(rhs.mDistTree)
805 , mDistAccessor(mDistTree)
806 , mIndexTree(rhs.mIndexTree)
807 , mIndexAccessor(mIndexTree)
808 , mIntersectionTree(rhs.mIntersectionTree)
809 , mIntersectionAccessor(mIntersectionTree)
814 template<
typename DistTreeT>
817 const tbb::blocked_range<size_t>& range)
const
819 typename StencilTreeT::LeafNodeType::ValueOnCIter iter;
821 for (
size_t n = range.begin(); n < range.end(); ++n) {
822 iter =
mLeafs.leaf(n).cbeginValueOn();
823 for (; iter; ++iter) {
824 evalVoxel(iter.getCoord());
829 template<
typename DistTreeT>
833 const DistValueT val = mDistAccessor.getValue(ijk),
zeroVal(0.0);
837 Vec3d dir = getClosestPointDir(ijk), n_dir;
842 for (
Int32 n = 0; n < 26; ++n) {
843 n_ijk = ijk + util::COORD_OFFSETS[n];
845 if (mIntersectionAccessor.isValueOn(n_ijk))
continue;
846 if (!mDistAccessor.probeValue(n_ijk, n_val))
continue;
849 n_dir = getClosestPointDir(n_ijk);
851 if (n_dir.dot(dir) > 0.0 ) {
853 mDistAccessor.setValue(ijk, -val);
859 template<
typename DistTreeT>
861 IntersectingVoxelSign<DistTreeT>::getClosestPointDir(
const Coord& ijk)
const
863 Vec3d voxelCenter(ijk[0], ijk[1], ijk[2]);
864 Vec4I prim = (*mPolygonList)[mIndexAccessor.getValue(ijk)];
867 Vec3d p0((*mPointList)[prim[0]]);
868 Vec3d p1((*mPointList)[prim[1]]);
869 Vec3d p2((*mPointList)[prim[2]]);
876 Vec3d p3((*mPointList)[prim[3]]);
887 Vec3d closestPoint = p0 * uv[0] +
889 p2 * (1.0 - uv(0) - uv(1));
891 Vec3d dir = (voxelCenter-closestPoint);
903 template<
typename DistTreeT>
911 typedef typename DistTreeT::template ValueConverter<Int32>::Type
IndexTreeT;
913 typedef typename DistTreeT::template ValueConverter<bool>::Type
StencilTreeT;
927 void operator()(
const tbb::blocked_range<size_t>&)
const;
932 DistTreeT& mDistTree;
943 template<
typename DistTreeT>
947 tbb::parallel_for(
mLeafs.getRange(), *
this);
948 mIntersectionTree.pruneInactive();
951 template<
typename DistTreeT>
955 (*this)(
mLeafs.getRange());
956 mIntersectionTree.pruneInactive();
959 template<
typename DistTreeT>
965 : mDistTree(distTree)
966 , mDistAccessor(mDistTree)
967 , mIndexTree(indexTree)
968 , mIndexAccessor(mIndexTree)
969 , mIntersectionTree(intersectionTree)
970 , mIntersectionAccessor(mIntersectionTree)
975 template<
typename DistTreeT>
978 : mDistTree(rhs.mDistTree)
979 , mDistAccessor(mDistTree)
980 , mIndexTree(rhs.mIndexTree)
981 , mIndexAccessor(mIndexTree)
982 , mIntersectionTree(rhs.mIntersectionTree)
983 , mIntersectionAccessor(mIntersectionTree)
988 template<
typename DistTreeT>
991 const tbb::blocked_range<size_t>& range)
const
995 DistValueT value, bg = mDistTree.getBackground();
997 typename StencilTreeT::LeafNodeType::ValueOnCIter iter;
999 for (
size_t n = range.begin(); n < range.end(); ++n) {
1001 typename StencilTreeT::LeafNodeType& leaf =
mLeafs.leaf(n);
1002 iter = leaf.cbeginValueOn();
1004 for (; iter; ++iter) {
1006 ijk = iter.getCoord();
1009 for (
Int32 m = 0; m < 26; ++m) {
1010 m_ijk = ijk + util::COORD_OFFSETS[m];
1011 if (mDistAccessor.probeValue(m_ijk, value)) {
1019 if (turnOff) leaf.setValueOff(ijk, bg);
1031 template<
typename DistTreeT>
1040 typedef typename DistTreeT::template ValueConverter<Int32>::Type
IndexTreeT;
1042 typedef typename DistTreeT::template ValueConverter<bool>::Type
StencilTreeT;
1055 void operator()(
const tbb::blocked_range<size_t>&)
const;
1060 DistTreeT& mDistTree;
1071 template<
typename DistTreeT>
1075 tbb::parallel_for(
mLeafs.getRange(), *
this);
1076 mDistTree.pruneInactive();
1077 mIndexTree.pruneInactive();
1080 template<
typename DistTreeT>
1084 (*this)(
mLeafs.getRange());
1085 mDistTree.pruneInactive();
1086 mIndexTree.pruneInactive();
1089 template<
typename DistTreeT>
1091 DistTreeT& distTree,
1095 : mDistTree(distTree)
1097 , mDistAccessor(mDistTree)
1098 , mIndexTree(indexTree)
1099 , mIndexAccessor(mIndexTree)
1100 , mIntersectionTree(intersectionTree)
1101 , mIntersectionAccessor(mIntersectionTree)
1105 template<
typename DistTreeT>
1108 : mDistTree(rhs.mDistTree)
1110 , mDistAccessor(mDistTree)
1111 , mIndexTree(rhs.mIndexTree)
1112 , mIndexAccessor(mIndexTree)
1113 , mIntersectionTree(rhs.mIntersectionTree)
1114 , mIntersectionAccessor(mIntersectionTree)
1118 template<
typename DistTreeT>
1121 const tbb::blocked_range<size_t>& range)
const
1128 const DistValueT distC = -0.86602540378443861;
1129 const DistValueT distBG = mDistTree.getBackground();
1130 const Int32 indexBG = mIntersectionTree.getBackground();
1132 typename DistTreeT::LeafNodeType::ValueOnCIter iter;
1133 for (
size_t n = range.begin(); n < range.end(); ++n) {
1135 typename DistTreeT::LeafNodeType& leaf =
mLeafs.leaf(n);
1136 iter = leaf.cbeginValueOn();
1137 for (; iter; ++iter) {
1139 value = iter.getValue();
1140 if(value > 0.0)
continue;
1142 ijk = iter.getCoord();
1143 if (mIntersectionAccessor.isValueOn(ijk))
continue;
1146 for (
Int32 m = 0; m < 18; ++m) {
1147 m_ijk = ijk + util::COORD_OFFSETS[m];
1149 if (mIntersectionAccessor.isValueOn(m_ijk)) {
1156 leaf.setValueOff(ijk, distBG);
1159 mIndexAccessor.setValueOff(ijk, indexBG);
1162 if (value > distC) leaf.setValue(ijk, distC);
1175 template<
typename DistTreeT>
1183 typedef typename DistTreeT::template ValueConverter<Int32>::Type
IndexTreeT;
1185 typedef typename DistTreeT::template ValueConverter<bool>::Type
StencilTreeT;
1190 ExpandNB(
const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList,
1201 void operator()(
const tbb::blocked_range<size_t>&)
const;
1206 double getDistToPrim(
const Coord& ijk,
const Int32 polyIdx)
const;
1208 std::vector<Vec3s>
const *
const mPointList;
1209 std::vector<Vec4I>
const *
const mPolygonList;
1211 DistTreeT& mDistTree;
1216 const DistValueT mExteriorBandWidth, mInteriorBandWidth, mVoxelSize;
1219 template<
typename DistTreeT>
1223 tbb::parallel_for(
mLeafs.getRange(), *
this);
1224 mMaskTree.pruneInactive();
1227 template<
typename DistTreeT>
1231 (*this)(
mLeafs.getRange());
1232 mMaskTree.pruneInactive();
1235 template<
typename DistTreeT>
1237 const std::vector<Vec3s>& pointList,
1238 const std::vector<Vec4I>& polygonList,
1239 DistTreeT& distTree,
1245 : mPointList(&pointList)
1246 , mPolygonList(&polygonList)
1247 , mDistTree(distTree)
1248 , mIndexTree(indexTree)
1249 , mMaskTree(maskTree)
1251 , mExteriorBandWidth(exteriorBandWidth)
1252 , mInteriorBandWidth(interiorBandWidth)
1253 , mVoxelSize(voxelSize)
1257 template<
typename DistTreeT>
1259 : mPointList(rhs.mPointList)
1260 , mPolygonList(rhs.mPolygonList)
1261 , mDistTree(rhs.mDistTree)
1262 , mIndexTree(rhs.mIndexTree)
1263 , mMaskTree(rhs.mMaskTree)
1265 , mExteriorBandWidth(rhs.mExteriorBandWidth)
1266 , mInteriorBandWidth(rhs.mInteriorBandWidth)
1267 , mVoxelSize(rhs.mVoxelSize)
1271 template<
typename DistTreeT>
1275 typedef typename DistTreeT::LeafNodeType DistLeafT;
1276 typedef typename IndexTreeT::LeafNodeType IndexLeafT;
1277 typedef typename StencilTreeT::LeafNodeType StencilLeafT;
1280 Int32 closestPrimIndex = 0;
1285 for (
size_t n = range.begin(); n < range.end(); ++n) {
1287 StencilLeafT& maskLeaf =
mLeafs.leaf(n);
1289 const Coord origin = maskLeaf.getOrigin();
1291 DistLeafT* distLeafPt = distAcc.
probeLeaf(origin);
1292 IndexLeafT* indexLeafPt = indexAcc.
probeLeaf(origin);
1294 if (distLeafPt != NULL && indexLeafPt != NULL) {
1296 DistLeafT& distLeaf = *distLeafPt;
1297 IndexLeafT& indexLeaf = *indexLeafPt;
1300 typename StencilLeafT::ValueOnIter iter = maskLeaf.beginValueOn();
1301 for (; iter; ++iter) {
1303 const Index pos = iter.pos();
1305 if (distLeaf.isValueOn(pos)) {
1311 getDist(iter.getCoord(), distAcc, indexAcc, closestPrimIndex);
1313 bool inside = distLeaf.getValue(pos) <
DistValueT(0.0);
1315 if (!inside && !(distance > mExteriorBandWidth)) {
1316 distLeaf.setValueOn(pos, distance);
1317 indexLeaf.setValueOn(pos, closestPrimIndex);
1318 }
else if (inside && !(distance > mInteriorBandWidth)) {
1319 distLeaf.setValueOn(pos, -distance);
1320 indexLeaf.setValueOn(pos, closestPrimIndex);
1327 maskLeaf.setValuesOff();
1332 template<
typename DistTreeT>
1335 IndexAccessorT& indexAcc,
Int32& primIndex)
const
1337 Vec3d voxelCenter(ijk[0], ijk[1], ijk[2]);
1342 for (
Int32 n = 0; n < 18; ++n) {
1343 n_ijk = ijk + util::COORD_OFFSETS[n];
1344 if (distAcc.probeValue(n_ijk, nDist)) {
1345 nDist = std::abs(nDist);
1348 primIndex = indexAcc.getValue(n_ijk);
1354 return getDistToPrim(ijk, primIndex);
1358 template<
typename DistTreeT>
1360 ExpandNB<DistTreeT>::getDistToPrim(
const Coord& ijk,
const Int32 polyIdx)
const
1362 Vec3d voxelCenter(ijk[0], ijk[1], ijk[2]);
1363 Vec4I verts = (*mPolygonList)[polyIdx];
1366 Vec3d p0((*mPointList)[verts[0]]);
1367 Vec3d p1((*mPointList)[verts[1]]);
1368 Vec3d p2((*mPointList)[verts[2]]);
1374 p1 =
Vec3d((*mPointList)[verts[3]]);
1377 if (secondDist < dist) dist = secondDist;
1380 return std::sqrt(dist) * double(mVoxelSize);
1395 template<
typename DistTreeT>
1398 typename DistTreeT::template ValueConverter<bool>::Type& intersectionTree)
1400 typedef typename DistTreeT::template ValueConverter<bool>::Type StencilTreeT;
1403 typedef typename DistTreeT::ValueType DistValueT;
1405 StencilAccessorT intrAccessor(intersectionTree);
1406 DistAccessorT distAccessor(distTree);
1408 std::deque<Coord> coordList;
1409 coordList.push_back(seed);
1412 while (!coordList.empty()) {
1413 ijk = coordList.back();
1414 coordList.pop_back();
1416 if (!distAccessor.isValueOn(ijk))
continue;
1418 DistValueT& dist =
const_cast<DistValueT&
>(distAccessor.getValue(ijk));
1419 if (!(dist < 0.0))
continue;
1422 for (
int n = 0; n < 6; ++n) {
1423 n_ijk = ijk + util::COORD_OFFSETS[n];
1425 if (!intrAccessor.isValueOn(n_ijk)) {
1426 if (distAccessor.isValueOn(n_ijk)) {
1427 if (distAccessor.getValue(n_ijk) < 0.0) {
1428 coordList.push_back(n_ijk);
1445 template<
typename DistTreeT,
typename InterruptT>
1448 typename DistTreeT::template ValueConverter<bool>::Type& intersectionTree,
1449 InterruptT *interrupter = NULL)
1451 typedef typename DistTreeT::template ValueConverter<bool>::Type StencilTreeT;
1454 typedef typename DistTreeT::ValueType DistValueT;
1456 StencilAccessorT intrAccessor(intersectionTree);
1457 DistAccessorT distAccessor(distTree);
1460 typename DistTreeT::LeafIter leafIter = distTree.beginLeaf();
1461 for (; leafIter; leafIter.next()) {
1463 if (interrupter && interrupter->wasInterrupted())
break;
1465 typename DistTreeT::LeafNodeType::ValueOnIter iter = leafIter->beginValueOn();
1466 for (; iter; iter.next()) {
1468 ijk = iter.getCoord();
1471 if (intrAccessor.isValueOn(ijk))
continue;
1473 if (iter.getValue() < 0.0) {
1474 for (
Int32 n = 0; n < 6; ++n) {
1475 n_ijk = ijk + util::COORD_OFFSETS[n];
1477 if (distAccessor.isValueOn(n_ijk) && distAccessor.getValue(n_ijk) > 0.0) {
1489 template<
typename ValueType>
1493 : mVoxelSize(voxelSize)
1494 , mUnsigned(unsignedDist)
1498 template <
typename LeafNodeType>
1505 typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
1506 for (; iter; ++iter) {
1507 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
1508 val = w[!mUnsigned && int(val < ValueType(0.0))] * std::sqrt(std::abs(val));
1513 ValueType mVoxelSize;
1514 const bool mUnsigned;
1518 template<
typename ValueType>
1522 : mExBandWidth(exBandWidth)
1523 , mInBandWidth(inBandWidth)
1527 template <
typename LeafNodeType>
1530 ValueType bgValues[2];
1531 bgValues[0] = mExBandWidth;
1532 bgValues[1] = -mInBandWidth;
1534 typename LeafNodeType::ValueOffIter iter = leaf.beginValueOff();
1536 for (; iter; ++iter) {
1537 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
1538 val = bgValues[int(val < ValueType(0.0))];
1543 ValueType mExBandWidth, mInBandWidth;
1547 template<
typename ValueType>
1550 TrimOp(ValueType exBandWidth, ValueType inBandWidth)
1551 : mExBandWidth(exBandWidth)
1552 , mInBandWidth(inBandWidth)
1556 template <
typename LeafNodeType>
1559 typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
1561 for (; iter; ++iter) {
1562 ValueType val = iter.getValue();
1563 if (val < -mInBandWidth) {
1564 iter.setValue(-mInBandWidth);
1566 }
else if (val > mExBandWidth) {
1567 iter.setValue(mExBandWidth);
1574 ValueType mExBandWidth, mInBandWidth;
1578 template<
typename ValueType>
1585 template <
typename LeafNodeType>
1588 typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
1589 for (; iter; ++iter) {
1590 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
1600 template<
typename Gr
idType,
typename ValueType>
1604 typedef typename Scheme::template ISStencil<GridType>::StencilType
Stencil;
1611 , mVoxelSize(voxelSize)
1618 template <
typename LeafNodeType>
1621 const ValueType dt = mCFL * mVoxelSize, one(1.0), invDx = one / mVoxelSize;
1626 typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
1627 for (; iter; ++iter) {
1628 stencil.moveTo(iter);
1630 const ValueType normSqGradPhi =
1633 const ValueType phi0 = stencil.getValue();
1634 const ValueType diff =
math::Sqrt(normSqGradPhi) * invDx - one;
1637 buffer.setValue(iter.pos(), phi0 - dt * S * diff);
1644 ValueType mVoxelSize, mCFL;
1648 template<
typename TreeType,
typename ValueType>
1656 template <
typename LeafNodeType>
1661 typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
1662 for (; iter; ++iter) {
1663 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
1664 val =
std::min(val, buffer.getValue(iter.pos()));
1673 template<
typename TreeType,
typename ValueType>
1681 , mBufferIndex(bufferIndex)
1685 template <
typename LeafNodeType>
1690 typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
1691 for (; iter; ++iter) {
1692 leaf.setValueOnly(iter.pos(), buffer.getValue(iter.pos()));
1698 const size_t mBufferIndex;
1709 template<
typename DistGr
idT,
typename InterruptT>
1711 openvdb::math::Transform::Ptr& transform,
int conversionFlags,
1712 InterruptT *interrupter,
int signSweeps)
1713 : mTransform(transform)
1714 , mConversionFlags(conversionFlags)
1715 , mSignSweeps(signSweeps)
1716 , mInterrupter(interrupter)
1719 mSignSweeps =
std::min(mSignSweeps, 1);
1723 template<
typename DistGr
idT,
typename InterruptT>
1729 mIntersectingVoxelsGrid = StencilGridT::create(
false);
1733 template<
typename DistGr
idT,
typename InterruptT>
1736 const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList,
1741 const DistValueT vs = mTransform->voxelSize()[0];
1742 doConvert(pointList, polygonList, vs * exBandWidth, vs * inBandWidth);
1747 template<
typename DistGr
idT,
typename InterruptT>
1750 const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList,
1754 const DistValueT vs = mTransform->voxelSize()[0];
1755 doConvert(pointList, polygonList, vs * exBandWidth, 0.0,
true);
1759 template<
typename DistGr
idT,
typename InterruptT>
1762 const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList,
1765 doConvert(pointList, polygonList, exBandWidth, inBandWidth);
1769 template<
typename DistGr
idT,
typename InterruptT>
1772 const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList,
1773 DistValueT exBandWidth, DistValueT inBandWidth,
bool unsignedDistField)
1775 mDistGrid->setTransform(mTransform);
1776 mIndexGrid->setTransform(mTransform);
1778 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1783 voxelizer(pointList, polygonList, mInterrupter);
1785 voxelizer.runParallel();
1787 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1789 mDistGrid->tree().merge(voxelizer.sqrDistTree());
1790 mIndexGrid->tree().merge(voxelizer.primIndexTree());
1791 mIntersectingVoxelsGrid->tree().merge(voxelizer.intersectionTree());
1794 if (!unsignedDistField) {
1799 internal::ContourTracer<DistTreeT, InterruptT> trace(
1800 mDistGrid->tree(), mIntersectingVoxelsGrid->tree(), mInterrupter);
1802 for (
int i = 0; i < mSignSweeps; ++i) {
1804 if (mInterrupter && mInterrupter->wasInterrupted())
break;
1805 trace.runParallel();
1807 if (mInterrupter && mInterrupter->wasInterrupted())
break;
1810 internal::propagateSign<DistTreeT, InterruptT>
1811 (mDistGrid->tree(), mIntersectingVoxelsGrid->tree(), mInterrupter);
1815 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1818 tree::LeafManager<StencilTreeT> leafs(mIntersectingVoxelsGrid->tree());
1821 internal::IntersectingVoxelSign<DistTreeT> sign(pointList, polygonList,
1822 mDistGrid->tree(), mIndexGrid->tree(), mIntersectingVoxelsGrid->tree(), leafs);
1826 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1830 internal::IntersectingVoxelCleaner<DistTreeT> cleaner(mDistGrid->tree(),
1831 mIndexGrid->tree(), mIntersectingVoxelsGrid->tree(), leafs);
1833 cleaner.runParallel();
1836 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1842 tree::LeafManager<DistTreeT> leafs(mDistGrid->tree());
1844 internal::ShellVoxelCleaner<DistTreeT> cleaner(mDistGrid->tree(),
1845 leafs, mIndexGrid->tree(), mIntersectingVoxelsGrid->tree());
1847 cleaner.runParallel();
1850 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1853 inBandWidth = DistValueT(0.0);
1856 if (mDistGrid->activeVoxelCount() == 0)
return;
1858 const DistValueT voxelSize(mTransform->voxelSize()[0]);
1862 typedef internal::SqrtAndScaleOp<DistValueT> XForm;
1863 tree::LeafManager<DistTreeT> leafs(mDistGrid->tree());
1864 XForm op(voxelSize, unsignedDistField);
1865 LeafTransformer<DistTreeT, XForm> transform(leafs, op);
1867 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1869 transform.runParallel();
1872 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1874 if (!unsignedDistField) {
1876 mDistGrid->tree().signedFloodFill();
1878 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1882 tree::LeafManager<DistTreeT> leafs(mDistGrid->tree());
1884 typedef internal::VoxelSignOp<DistValueT> SignXForm;
1885 SignXForm op(exBandWidth, inBandWidth);
1887 LeafTransformer<DistTreeT, SignXForm> transform(leafs, op);
1888 transform.runParallel();
1890 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1892 DistValueT bgValues[2];
1893 bgValues[0] = exBandWidth;
1894 bgValues[1] = -inBandWidth;
1896 typename DistTreeT::ValueAllIter tileIt(mDistGrid->tree());
1897 tileIt.setMaxDepth(DistTreeT::ValueAllIter::LEAF_DEPTH - 1);
1899 for ( ; tileIt; ++tileIt) {
1900 DistValueT& val =
const_cast<DistValueT&
>(tileIt.getValue());
1901 val = bgValues[int(val < DistValueT(0.0))];
1904 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1907 typename DistTreeT::Ptr newTree(
new DistTreeT(exBandWidth));
1908 newTree->merge(mDistGrid->tree());
1909 mDistGrid->setTree(newTree);
1915 typedef internal::OffsetOp<DistValueT> OffsetOp;
1916 typedef internal::RenormOp<DistGridT, DistValueT> RenormOp;
1917 typedef internal::MinOp<DistTreeT, DistValueT> MinOp;
1918 typedef internal::MergeBufferOp<DistTreeT, DistValueT> MergeBufferOp;
1920 tree::LeafManager<DistTreeT> leafs(mDistGrid->tree(), 1);
1922 const DistValueT offset = 0.8 * voxelSize;
1924 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1926 OffsetOp offsetOp(-offset);
1927 LeafTransformer<DistTreeT, OffsetOp> offsetXform(leafs, offsetOp);
1928 offsetXform.runParallel();
1930 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1932 RenormOp renormOp(*mDistGrid, leafs, voxelSize);
1933 LeafTransformer<DistTreeT, RenormOp> renormXform(leafs, renormOp);
1934 renormXform.runParallel();
1937 LeafTransformer<DistTreeT, MinOp> minXform(leafs, minOp);
1938 minXform.runParallel();
1940 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1942 offsetOp.resetOffset(offset);
1943 offsetXform.runParallel();
1946 mIntersectingVoxelsGrid->clear();
1949 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1952 if (inBandWidth > voxelSize || exBandWidth > voxelSize) {
1955 StencilTreeT maskTree(
false);
1956 tree::ValueAccessor<StencilTreeT> acc(maskTree);
1957 maskTree.topologyUnion(mDistGrid->tree());
1961 typedef typename DistTreeT::LeafNodeType DistLeafType;
1963 std::vector<DistLeafType*> distLeafs;
1964 distLeafs.reserve(mDistGrid->tree().leafCount());
1966 typename DistTreeT::LeafIter iter = mDistGrid->tree().beginLeaf();
1967 for ( ; iter; ++iter) distLeafs.push_back(iter.getLeaf());
1969 tree::ValueAccessor<DistTreeT> distAcc(mDistGrid->tree());
1971 DistValueT leafSize = DistValueT(DistLeafType::DIM - 1) * voxelSize;
1973 const double inLeafsRatio = double(inBandWidth) / double(leafSize);
1975 if (
double(inLeafs) > (inLeafsRatio + 1.0)) {
1976 inLeafs = size_t(std::ceil(inLeafsRatio)) + 1;
1978 size_t exLeafs = size_t(std::ceil(exBandWidth / leafSize)) + 1;
1979 size_t numLeafs =
std::max(inLeafs, exLeafs);
1981 for (
size_t i = 0; i < numLeafs; ++i) {
1983 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1985 std::vector<DistLeafType*> newDistLeafs;
1986 newDistLeafs.reserve(2 * distLeafs.size());
1988 for (
size_t n = 0, N = distLeafs.size(); n < N; ++n) {
1990 Coord ijk = distLeafs[n]->getOrigin();
1992 const bool inside = distLeafs[n]->getValue(ijk) < DistValueT(0.0);
1994 if (inside && !(i < inLeafs))
continue;
1995 else if (!inside && !(i < exLeafs))
continue;
1998 if (distAcc.probeLeaf(ijk) == NULL) {
1999 newDistLeafs.push_back(distAcc.touchLeaf(ijk));
2004 if (distAcc.probeLeaf(ijk) == NULL) {
2005 newDistLeafs.push_back(distAcc.touchLeaf(ijk));
2010 if (distAcc.probeLeaf(ijk) == NULL) {
2011 newDistLeafs.push_back(distAcc.touchLeaf(ijk));
2015 ijk[0] += DistLeafType::DIM;
2016 if (distAcc.probeLeaf(ijk) == NULL) {
2017 newDistLeafs.push_back(distAcc.touchLeaf(ijk));
2020 ijk[0] -= DistLeafType::DIM;
2021 ijk[1] += DistLeafType::DIM;
2022 if (distAcc.probeLeaf(ijk) == NULL) {
2023 newDistLeafs.push_back(distAcc.touchLeaf(ijk));
2026 ijk[1] -= DistLeafType::DIM;
2027 ijk[2] += DistLeafType::DIM;
2028 if (distAcc.probeLeaf(ijk) == NULL) {
2029 newDistLeafs.push_back(distAcc.touchLeaf(ijk));
2033 if (newDistLeafs.empty())
break;
2034 distLeafs.swap(newDistLeafs);
2038 if (mInterrupter && mInterrupter->wasInterrupted())
return;
2040 mIndexGrid->tree().topologyUnion(mDistGrid->tree());
2042 while (maskTree.activeVoxelCount() > 0) {
2044 if (mInterrupter && mInterrupter->wasInterrupted())
break;
2047 tree::LeafManager<StencilTreeT> leafs(maskTree);
2049 internal::ExpandNB<DistTreeT> expand(pointList, polygonList, mDistGrid->tree(),
2050 mIndexGrid->tree(), maskTree, leafs, exBandWidth, inBandWidth, voxelSize);
2052 expand.runParallel();
2060 if ((inBandWidth < (voxelSize*3)) || (exBandWidth < (voxelSize*3))) {
2066 tree::LeafManager<DistTreeT> leafs(mDistGrid->tree());
2068 typedef internal::TrimOp<DistValueT> TrimOp;
2070 TrimOp op(exBandWidth, inBandWidth);
2071 LeafTransformer<DistTreeT, TrimOp> transform(leafs, op);
2072 transform.runParallel();
2075 if (mInterrupter && mInterrupter->wasInterrupted())
return;
2077 tree::LevelSetPrune<DistValueT> prune;
2078 mDistGrid->tree().pruneOp(prune);
2085 #endif // OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED