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/Morphology.h>
39 #include <openvdb/util/NullInterrupter.h>
40 #include <openvdb/util/Util.h>
42 #include <tbb/blocked_range.h>
43 #include <tbb/parallel_for.h>
44 #include <tbb/parallel_reduce.h>
76 template<
typename Gr
idType>
77 inline typename GridType::Ptr
79 const openvdb::math::Transform& xform,
80 const std::vector<Vec3s>& points,
81 const std::vector<Vec3I>& triangles,
100 template<
typename Gr
idType>
101 inline typename GridType::Ptr
103 const openvdb::math::Transform& xform,
104 const std::vector<Vec3s>& points,
105 const std::vector<Vec4I>& quads,
125 template<
typename Gr
idType>
126 inline typename GridType::Ptr
128 const openvdb::math::Transform& xform,
129 const std::vector<Vec3s>& points,
130 const std::vector<Vec3I>& triangles,
131 const std::vector<Vec4I>& quads,
153 template<
typename Gr
idType>
154 inline typename GridType::Ptr
156 const openvdb::math::Transform& xform,
157 const std::vector<Vec3s>& points,
158 const std::vector<Vec3I>& triangles,
159 const std::vector<Vec4I>& quads,
178 template<
typename Gr
idType>
179 inline typename GridType::Ptr
181 const openvdb::math::Transform& xform,
182 const std::vector<Vec3s>& points,
183 const std::vector<Vec3I>& triangles,
184 const std::vector<Vec4I>& quads,
196 template<
typename FloatGr
idT,
typename InterruptT = util::NullInterrupter>
202 typedef typename FloatTreeT::template ValueConverter<Int32>::Type
IntTreeT;
204 typedef typename FloatTreeT::template ValueConverter<bool>::Type
BoolTreeT;
207 MeshToVolume(openvdb::math::Transform::Ptr&,
int conversionFlags = 0,
208 InterruptT *interrupter = NULL,
int signSweeps = 1);
221 void convertToLevelSet(
222 const std::vector<Vec3s>& pointList,
223 const std::vector<Vec4I>& polygonList,
235 void convertToUnsignedDistanceField(
const std::vector<Vec3s>& pointList,
236 const std::vector<Vec4I>& polygonList,
FloatValueT exBandWidth);
241 typename FloatGridT::Ptr
distGridPtr()
const {
return mDistGrid; }
251 void doConvert(
const std::vector<Vec3s>&,
const std::vector<Vec4I>&,
252 FloatValueT exBandWidth, FloatValueT inBandWidth,
bool unsignedDistField =
false);
255 return mInterrupter && mInterrupter->wasInterrupted(percent);
258 openvdb::math::Transform::Ptr mTransform;
259 int mConversionFlags, mSignSweeps;
261 typename FloatGridT::Ptr mDistGrid;
262 typename IntGridT::Ptr mIndexGrid;
263 typename BoolGridT::Ptr mIntersectingVoxelsGrid;
265 InterruptT *mInterrupter;
282 : mXDist(dist), mYDist(dist), mZDist(dist)
325 void convert(
const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList);
331 std::vector<Vec3d>& points, std::vector<Index32>& primitives);
356 PointTransform(
const std::vector<Vec3s>& pointsIn, std::vector<Vec3s>& pointsOut,
358 : mPointsIn(pointsIn)
359 , mPointsOut(&pointsOut)
366 tbb::parallel_for(tbb::blocked_range<size_t>(0, mPointsOut->size()), *
this);
371 (*this)(tbb::blocked_range<size_t>(0, mPointsOut->size()));
374 inline void operator()(
const tbb::blocked_range<size_t>& range)
const
376 for (
size_t n = range.begin(); n < range.end(); ++n) {
377 (*mPointsOut)[n] = mXform.worldToIndex(mPointsIn[n]);
382 const std::vector<Vec3s>& mPointsIn;
383 std::vector<Vec3s> *
const mPointsOut;
388 template<
typename ValueType>
391 static ValueType
epsilon() {
return ValueType(1e-7); }
396 template<
typename FloatTreeT,
typename IntTreeT>
398 combine(FloatTreeT& lhsDist, IntTreeT& lhsIndex, FloatTreeT& rhsDist, IntTreeT& rhsIndex)
400 typedef typename FloatTreeT::ValueType FloatValueT;
404 typename FloatTreeT::LeafCIter iter = rhsDist.cbeginLeaf();
406 FloatValueT rhsValue;
409 for ( ; iter; ++iter) {
410 typename FloatTreeT::LeafNodeType::ValueOnCIter it = iter->cbeginValueOn();
415 rhsValue = it.getValue();
416 FloatValueT& lhsValue =
const_cast<FloatValueT&
>(lhsDistAccessor.
getValue(ijk));
418 if (-rhsValue < std::abs(lhsValue)) {
437 template<
typename FloatTreeT,
typename InterruptT = util::NullInterrupter>
444 typedef typename FloatTreeT::template ValueConverter<Int32>::Type
IntTreeT;
447 typedef typename FloatTreeT::template ValueConverter<bool>::Type
BoolTreeT;
453 const std::vector<Vec4I>& polygonList, InterruptT *interrupter = NULL);
457 void run(
bool threaded =
true);
460 void operator() (
const tbb::blocked_range<size_t> &range);
470 bool wasInterrupted()
const {
return mInterrupter && mInterrupter->wasInterrupted(); }
472 bool evalVoxel(
const Coord& ijk,
const Int32 polyIdx);
474 const std::vector<Vec3s>& mPointList;
475 const std::vector<Vec4I>& mPolygonList;
477 FloatTreeT mSqrDistTree;
478 FloatAccessorT mSqrDistAccessor;
480 IntTreeT mPrimIndexTree;
481 IntAccessorT mPrimIndexAccessor;
483 BoolTreeT mIntersectionTree;
484 BoolAccessorT mIntersectionAccessor;
487 IntTreeT mLastPrimTree;
488 IntAccessorT mLastPrimAccessor;
490 InterruptT *mInterrupter;
493 struct Primitive {
Vec3d a, b, c, d;
Int32 index; };
495 template<
bool IsQuad>
496 bool evalPrimitive(
const Coord&,
const Primitive&);
498 template<
bool IsQuad>
499 void voxelize(
const Primitive&);
503 template<
typename FloatTreeT,
typename InterruptT>
508 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPolygonList.size()), *
this);
510 (*this)(tbb::blocked_range<size_t>(0, mPolygonList.size()));
514 template<
typename FloatTreeT,
typename InterruptT>
516 const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList,
517 InterruptT *interrupter)
518 : mPointList(pointList)
519 , mPolygonList(polygonList)
521 , mSqrDistAccessor(mSqrDistTree)
523 , mPrimIndexAccessor(mPrimIndexTree)
524 , mIntersectionTree(false)
525 , mIntersectionAccessor(mIntersectionTree)
527 , mLastPrimAccessor(mLastPrimTree)
528 , mInterrupter(interrupter)
532 template<
typename FloatTreeT,
typename InterruptT>
535 : mPointList(rhs.mPointList)
536 , mPolygonList(rhs.mPolygonList)
538 , mSqrDistAccessor(mSqrDistTree)
540 , mPrimIndexAccessor(mPrimIndexTree)
541 , mIntersectionTree(false)
542 , mIntersectionAccessor(mIntersectionTree)
544 , mLastPrimAccessor(mLastPrimTree)
545 , mInterrupter(rhs.mInterrupter)
550 template<
typename FloatTreeT,
typename InterruptT>
556 for (
size_t n = range.begin(); n < range.end(); ++n) {
558 if (mInterrupter && mInterrupter->wasInterrupted()) {
559 tbb::task::self().cancel_group_execution();
563 const Vec4I& verts = mPolygonList[n];
565 prim.index =
Int32(n);
566 prim.a =
Vec3d(mPointList[verts[0]]);
567 prim.b =
Vec3d(mPointList[verts[1]]);
568 prim.c =
Vec3d(mPointList[verts[2]]);
571 prim.d =
Vec3d(mPointList[verts[3]]);
572 voxelize<true>(prim);
574 voxelize<false>(prim);
580 template<
typename FloatTreeT,
typename InterruptT>
581 template<
bool IsQuad>
585 std::deque<Coord> coordList;
589 coordList.push_back(ijk);
591 evalPrimitive<IsQuad>(ijk, prim);
593 while (!coordList.empty()) {
596 ijk = coordList.back();
597 coordList.pop_back();
599 mIntersectionAccessor.setActiveState(ijk,
true);
601 for (
Int32 i = 0; i < 26; ++i) {
604 if (prim.index != mLastPrimAccessor.getValue(nijk)) {
605 mLastPrimAccessor.setValue(nijk, prim.index);
606 if(evalPrimitive<IsQuad>(nijk, prim)) coordList.push_back(nijk);
613 template<
typename FloatTreeT,
typename InterruptT>
614 template<
bool IsQuad>
616 MeshVoxelizer<FloatTreeT, InterruptT>::evalPrimitive(
const Coord& ijk,
const Primitive& prim)
618 Vec3d uvw, voxelCenter(ijk[0], ijk[1], ijk[2]);
621 double dist = (voxelCenter -
626 double secondDist = (voxelCenter -
629 if (secondDist < dist) dist = secondDist;
632 FloatValueT oldDist = mSqrDistAccessor.getValue(ijk);
636 if (std::abs(oldDist) > dist) {
637 mSqrDistAccessor.setValue(ijk, -dist);
638 mPrimIndexAccessor.setValue(ijk, prim.index);
641 return (dist < 0.86602540378443861);
645 template<
typename FloatTreeT,
typename InterruptT>
649 typedef typename FloatTreeT::RootNodeType FloatRootNodeT;
651 typedef typename boost::mpl::at<FloatNodeTVec, boost::mpl::int_<1> >::type FloatInternalNodeT;
653 typedef typename IntTreeT::RootNodeType IntRootNodeT;
655 typedef typename boost::mpl::at<IntNodeTVec, boost::mpl::int_<1> >::type IntInternalNodeT;
660 rhs.mSqrDistTree.clearAllAccessors();
661 rhs.mPrimIndexTree.clearAllAccessors();
663 typename FloatTreeT::LeafIter leafIt = rhs.mSqrDistTree.beginLeaf();
664 for ( ; leafIt; ++leafIt) {
666 ijk = leafIt->getOrigin();
667 FloatLeafT* lhsDistLeafPt = mSqrDistAccessor.probeLeaf(ijk);
669 if (!lhsDistLeafPt) {
676 mSqrDistAccessor.addLeaf(rhs.mSqrDistAccessor.probeLeaf(ijk));
677 FloatInternalNodeT* floatNode =
678 rhs.mSqrDistAccessor.template getNode<FloatInternalNodeT>();
679 floatNode->template stealNode<FloatLeafT>(ijk,
FloatValueT(0.0),
false);
680 rhs.mSqrDistAccessor.clear();
682 mPrimIndexAccessor.addLeaf(rhs.mPrimIndexAccessor.probeLeaf(ijk));
683 IntInternalNodeT* intNode =
684 rhs.mPrimIndexAccessor.template getNode<IntInternalNodeT>();
685 intNode->template stealNode<IntLeafT>(ijk, 0,
false);
686 rhs.mPrimIndexAccessor.clear();
690 IntLeafT& lhsIdxLeaf = *mPrimIndexAccessor.probeLeaf(ijk);
693 IntLeafT& rhsIdxLeaf = *rhs.mPrimIndexAccessor.probeLeaf(ijk);
695 typename FloatLeafT::ValueOnCIter it = rhsDistLeaf.cbeginValueOn();
698 const FloatValueT& rhsValue = rhsDistLeaf.getValue(offset);
699 if (!lhsDistLeaf.isValueOn(offset)) {
700 lhsDistLeaf.setValueOn(offset, rhsValue);
701 lhsIdxLeaf.setValueOn(offset, rhsIdxLeaf.getValue(offset));
702 }
else if (rhsValue > lhsDistLeaf.getValue(offset)) {
703 lhsDistLeaf.setValueOnly(offset, rhsValue);
704 lhsIdxLeaf.setValueOnly(offset, rhsIdxLeaf.getValue(offset));
710 mIntersectionTree.merge(rhs.mIntersectionTree);
720 template<
typename FloatTreeT,
typename InterruptT = util::NullInterrupter>
726 typedef typename FloatTreeT::template ValueConverter<bool>::Type
BoolTreeT;
736 void operator()(
const tbb::blocked_range<int> &range)
const;
741 int sparseScan(
int slice)
const;
743 FloatTreeT& mDistTree;
752 std::vector<Index> mStepSize;
754 InterruptT *mInterrupter;
758 template<
typename FloatTreeT,
typename InterruptT>
762 tbb::parallel_for(tbb::blocked_range<int>(mBBox.min()[0], mBBox.max()[0]+1), *
this);
766 template<
typename FloatTreeT,
typename InterruptT>
770 (*this)(tbb::blocked_range<int>(mBBox.min()[0], mBBox.max()[0]+1));
774 template<
typename FloatTreeT,
typename InterruptT>
776 FloatTreeT& distTree,
const BoolTreeT& intersectionTree, InterruptT *interrupter)
777 : mDistTree(distTree)
778 , mDistAccessor(mDistTree)
779 , mIntersectionTree(intersectionTree)
780 , mIntersectionAccessor(mIntersectionTree)
783 , mInterrupter(interrupter)
786 std::vector<Index> dims;
787 mDistTree.getNodeLog2Dims(dims);
789 mStepSize.resize(dims.size()+1, 1);
791 for (
int idx = static_cast<int>(dims.size()) - 1; idx > -1; --idx) {
792 exponent += dims[idx];
793 mStepSize[idx] = 1 << exponent;
796 mDistTree.evalLeafBoundingBox(mBBox);
799 const int tileDim = mStepSize[0];
801 for (
size_t i = 0; i < 3; ++i) {
804 double diff = std::abs(
double(mBBox.
min()[i])) / double(tileDim);
806 if (mBBox.
min()[i] <= tileDim) {
807 n = int(std::ceil(diff));
808 mBBox.
min()[i] = - n * tileDim;
810 n = int(std::floor(diff));
811 mBBox.
min()[i] = n * tileDim;
814 n = int(std::ceil(std::abs(
double(mBBox.
max()[i] - mBBox.
min()[i])) / double(tileDim)));
815 mBBox.
max()[i] = mBBox.
min()[i] + n * tileDim;
820 template<
typename FloatTreeT,
typename InterruptT>
823 : mDistTree(rhs.mDistTree)
824 , mDistAccessor(mDistTree)
825 , mIntersectionTree(rhs.mIntersectionTree)
826 , mIntersectionAccessor(mIntersectionTree)
828 , mStepSize(rhs.mStepSize)
829 , mInterrupter(rhs.mInterrupter)
834 template<
typename FloatTreeT,
typename InterruptT>
840 for (
int n = range.begin(); n < range.end(); n += iStep) {
842 if (mInterrupter && mInterrupter->wasInterrupted()) {
843 tbb::task::self().cancel_group_execution();
847 iStep = sparseScan(n);
852 template<
typename FloatTreeT,
typename InterruptT>
856 bool lastVoxelWasOut =
true;
859 Coord ijk(slice, mBBox.min()[1], mBBox.min()[2]);
860 Coord step(mStepSize[mDistAccessor.getValueDepth(ijk) + 1]);
863 for (ijk[1] = mBBox.min()[1]; ijk[1] <= mBBox.max()[1]; ijk[1] += step[1]) {
865 if (mInterrupter && mInterrupter->wasInterrupted()) {
869 step[1] = mStepSize[mDistAccessor.getValueDepth(ijk) + 1];
870 step[0] =
std::min(step[0], step[1]);
872 for (ijk[2] = mBBox.min()[2]; ijk[2] <= mBBox.max()[2]; ijk[2] += step[2]) {
874 step[2] = mStepSize[mDistAccessor.getValueDepth(ijk) + 1];
875 step[1] =
std::min(step[1], step[2]);
876 step[0] =
std::min(step[0], step[2]);
879 if (mDistAccessor.isValueOn(ijk)) {
882 if (mIntersectionAccessor.isValueOn(ijk)) {
884 lastVoxelWasOut =
false;
887 }
else if (lastVoxelWasOut) {
889 FloatValueT& val =
const_cast<FloatValueT&
>(mDistAccessor.getValue(ijk));
895 for (
Int32 n = 3; n < 6; n += 2) {
896 n_ijk = ijk + util::COORD_OFFSETS[n];
898 if (mDistAccessor.probeValue(n_ijk, val) && val > 0) {
899 lastVoxelWasOut =
true;
904 if (lastVoxelWasOut) {
906 FloatValueT& v =
const_cast<FloatValueT&
>(mDistAccessor.getValue(ijk));
909 const int tmp_k = ijk[2];
912 for (--ijk[2]; ijk[2] >= last_k; --ijk[2]) {
913 if (mIntersectionAccessor.isValueOn(ijk))
break;
915 const_cast<FloatValueT&
>(mDistAccessor.getValue(ijk));
916 if (vb < FloatValueT(0.0)) vb = -vb;
939 template<
typename FloatTreeT,
typename InterruptT = util::NullInterrupter>
947 typedef typename FloatTreeT::template ValueConverter<bool>::Type
BoolTreeT;
954 InterruptT *interrupter = NULL);
962 void operator() (
const tbb::blocked_range<size_t> &range);
970 bool wasInterrupted()
const {
return mInterrupter && mInterrupter->wasInterrupted(); }
973 const FloatTreeT& mDistTree;
978 InterruptT *mInterrupter;
982 template<
typename FloatTreeT,
typename InterruptT>
985 const BoolTreeT& intersectionTree, InterruptT *interrupter)
986 : mDistLeafs(distLeafs)
987 , mDistTree(distTree)
988 , mIntersectionTree(intersectionTree)
989 , mSignMaskTree(false)
990 , mInterrupter(interrupter)
995 template<
typename FloatTreeT,
typename InterruptT>
998 : mDistLeafs(rhs.mDistLeafs)
999 , mDistTree(rhs.mDistTree)
1000 , mIntersectionTree(rhs.mIntersectionTree)
1001 , mSignMaskTree(false)
1002 , mInterrupter(rhs.mInterrupter)
1007 template<
typename FloatTreeT,
typename InterruptT>
1011 tbb::parallel_reduce(mDistLeafs.getRange(), *
this);
1015 template<
typename FloatTreeT,
typename InterruptT>
1019 (*this)(mDistLeafs.getRange());
1023 template<
typename FloatTreeT,
typename InterruptT>
1036 const int extent = BoolLeafT::DIM - 1;
1038 for (
size_t n = range.begin(); n < range.end(); ++n) {
1040 const FloatLeafT& distLeaf = mDistLeafs.leaf(n);
1042 minCoord = distLeaf.getOrigin();
1043 maxCoord[0] = minCoord[0] + extent;
1044 maxCoord[1] = minCoord[1] + extent;
1045 maxCoord[2] = minCoord[2] + extent;
1051 bool addLeaf =
false;
1055 typename FloatLeafT::ValueOnCIter it = distLeaf.cbeginValueOn();
1057 if (intersectionLeaf && intersectionLeaf->isValueOn(it.pos()))
continue;
1059 ijk = it.getCoord();
1061 for (
size_t i = 0; i < 6; ++i) {
1062 if (distLeaf.probeValue(ijk+util::COORD_OFFSETS[i], value) && value>0.0) {
1063 maskLeaf.setValueOn(ijk);
1069 for (
size_t i = 0; i < 6; ++i) {
1070 if (distAcc.
probeValue(ijk+util::COORD_OFFSETS[i], value) && value>0.0) {
1071 maskLeaf.setValueOn(ijk);
1080 if (addLeaf) maskAcc.
addLeaf(maskLeafPt);
1081 else delete maskLeafPt;
1086 template<
typename FloatTreeT,
typename InterruptT>
1090 mSignMaskTree.merge(rhs.mSignMaskTree);
1098 template<
typename FloatTreeT,
typename InterruptT = util::NullInterrupter>
1105 typedef typename FloatTreeT::template ValueConverter<bool>::Type
BoolTreeT;
1119 void operator() (
const tbb::blocked_range<size_t> &range);
1127 bool wasInterrupted()
const {
return mInterrupter && mInterrupter->wasInterrupted(); }
1130 FloatTreeT& mDistTree;
1134 InterruptT *mInterrupter;
1138 template<
typename FloatTreeT,
typename InterruptT>
1140 FloatTreeT& distTree,
const BoolTreeT& intersectionTree, InterruptT *interrupter)
1141 : mOldSignMaskLeafs(signMaskLeafs)
1142 , mDistTree(distTree)
1143 , mIntersectionTree(intersectionTree)
1144 , mSignMaskTree(false)
1145 , mInterrupter(interrupter)
1150 template<
typename FloatTreeT,
typename InterruptT>
1153 : mOldSignMaskLeafs(rhs.mOldSignMaskLeafs)
1154 , mDistTree(rhs.mDistTree)
1155 , mIntersectionTree(rhs.mIntersectionTree)
1156 , mSignMaskTree(false)
1157 , mInterrupter(rhs.mInterrupter)
1162 template<
typename FloatTreeT,
typename InterruptT>
1166 tbb::parallel_reduce(mOldSignMaskLeafs.getRange(), *
this);
1170 template<
typename FloatTreeT,
typename InterruptT>
1174 (*this)(mOldSignMaskLeafs.getRange());
1178 template<
typename FloatTreeT,
typename InterruptT>
1186 std::deque<Coord> coordList;
1193 const int extent = BoolLeafT::DIM - 1;
1195 for (
size_t n = range.begin(); n < range.end(); ++n) {
1196 BoolLeafT& oldMaskLeaf = mOldSignMaskLeafs.leaf(n);
1198 minCoord = oldMaskLeaf.getOrigin();
1199 maxCoord[0] = minCoord[0] + extent;
1200 maxCoord[1] = minCoord[1] + extent;
1201 maxCoord[2] = minCoord[2] + extent;
1206 typename BoolLeafT::ValueOnCIter it = oldMaskLeaf.cbeginValueOn();
1208 coordList.push_back(it.getCoord());
1210 while (!coordList.empty()) {
1212 ijk = coordList.back();
1213 coordList.pop_back();
1219 for (
size_t i = 0; i < 6; ++i) {
1220 nijk = ijk + util::COORD_OFFSETS[i];
1222 if (intersectionLeaf && intersectionLeaf->isValueOn(nijk))
continue;
1224 if (distLeaf.probeValue(nijk, value) && value < 0.0) {
1225 coordList.push_back(nijk);
1230 distAcc.
probeValue(nijk, value) && value < 0.0) {
1242 template<
typename FloatTreeT,
typename InterruptT>
1246 mSignMaskTree.merge(rhs.mSignMaskTree);
1257 template<
typename FloatTreeT>
1263 typedef typename FloatTreeT::template ValueConverter<Int32>::Type
IntTreeT;
1265 typedef typename FloatTreeT::template ValueConverter<bool>::Type
BoolTreeT;
1270 const std::vector<Vec3s>& pointList,
1271 const std::vector<Vec4I>& polygonList,
1272 FloatTreeT& distTree,
1283 void operator()(
const tbb::blocked_range<size_t>&)
const;
1290 std::vector<Vec3s>
const *
const mPointList;
1291 std::vector<Vec4I>
const *
const mPolygonList;
1293 FloatTreeT& mDistTree;
1301 template<
typename FloatTreeT>
1305 tbb::parallel_for(mLeafs.getRange(), *
this);
1309 template<
typename FloatTreeT>
1313 (*this)(mLeafs.getRange());
1317 template<
typename FloatTreeT>
1319 const std::vector<Vec3s>& pointList,
1320 const std::vector<Vec4I>& polygonList,
1321 FloatTreeT& distTree,
1325 : mPointList(&pointList)
1326 , mPolygonList(&polygonList)
1327 , mDistTree(distTree)
1328 , mIndexTree(indexTree)
1329 , mIntersectionTree(intersectionTree)
1335 template<
typename FloatTreeT>
1338 : mPointList(rhs.mPointList)
1339 , mPolygonList(rhs.mPolygonList)
1340 , mDistTree(rhs.mDistTree)
1341 , mIndexTree(rhs.mIndexTree)
1342 , mIntersectionTree(rhs.mIntersectionTree)
1343 , mLeafs(rhs.mLeafs)
1348 template<
typename FloatTreeT>
1351 const tbb::blocked_range<size_t>& range)
const
1362 typename BoolTreeT::LeafNodeType::ValueOnCIter iter;
1363 for (
size_t n = range.begin(); n < range.end(); ++n) {
1364 iter = mLeafs.leaf(n).cbeginValueOn();
1365 for (; iter; ++iter) {
1367 ijk = iter.getCoord();
1373 center =
Vec3d(ijk[0], ijk[1], ijk[2]);
1375 for (
Int32 i = 0; i < 26; ++i) {
1376 nijk = ijk + util::COORD_OFFSETS[i];
1383 cpt = getClosestPoint(nijk, prim);
1385 dir1 = center -
cpt;
1388 dir2 =
Vec3d(nijk[0], nijk[1], nijk[2]) -
cpt;
1391 if (dir2.
dot(dir1) > 0.0) {
1402 template<
typename FloatTreeT>
1406 Vec3d voxelCenter(ijk[0], ijk[1], ijk[2]);
1409 const Vec3d a((*mPointList)[prim[0]]);
1410 const Vec3d b((*mPointList)[prim[1]]);
1411 const Vec3d c((*mPointList)[prim[2]]);
1419 Vec3d diff1 = voxelCenter - cpt1;
1421 const Vec3d d((*mPointList)[prim[3]]);
1424 Vec3d diff2 = voxelCenter - cpt2;
1441 template<
typename FloatTreeT>
1448 typedef typename FloatTreeT::template ValueConverter<Int32>::Type
IntTreeT;
1451 typedef typename FloatTreeT::template ValueConverter<bool>::Type
BoolTreeT;
1465 void operator()(
const tbb::blocked_range<size_t>&)
const;
1470 FloatTreeT& mDistTree;
1477 template<
typename FloatTreeT>
1481 tbb::parallel_for(mLeafs.getRange(), *
this);
1482 mIntersectionTree.pruneInactive();
1486 template<
typename FloatTreeT>
1490 (*this)(mLeafs.getRange());
1491 mIntersectionTree.pruneInactive();
1495 template<
typename FloatTreeT>
1497 FloatTreeT& distTree,
1501 : mDistTree(distTree)
1502 , mIndexTree(indexTree)
1503 , mIntersectionTree(intersectionTree)
1509 template<
typename FloatTreeT>
1512 : mDistTree(rhs.mDistTree)
1513 , mIndexTree(rhs.mIndexTree)
1514 , mIntersectionTree(rhs.mIntersectionTree)
1515 , mLeafs(rhs.mLeafs)
1520 template<
typename FloatTreeT>
1523 const tbb::blocked_range<size_t>& range)
const
1530 typename BoolLeafT::ValueOnCIter iter;
1536 for (
size_t n = range.begin(); n < range.end(); ++n) {
1540 ijk = maskLeaf.getOrigin();
1544 iter = maskLeaf.cbeginValueOn();
1545 for (; iter; ++iter) {
1547 offset = iter.pos();
1549 if(distLeaf->getValue(offset) > 0.1)
continue;
1551 ijk = iter.getCoord();
1553 for (
Int32 m = 0; m < 26; ++m) {
1554 m_ijk = ijk + util::COORD_OFFSETS[m];
1564 maskLeaf.setValueOff(offset);
1565 distLeaf->setValueOn(offset, -0.86602540378443861);
1579 template<
typename FloatTreeT>
1587 typedef typename FloatTreeT::template ValueConverter<Int32>::Type
IntTreeT;
1590 typedef typename FloatTreeT::template ValueConverter<bool>::Type
BoolTreeT;
1603 void operator()(
const tbb::blocked_range<size_t>&)
const;
1608 FloatTreeT& mDistTree;
1615 template<
typename FloatTreeT>
1619 tbb::parallel_for(mLeafs.getRange(), *
this);
1620 mDistTree.pruneInactive();
1621 mIndexTree.pruneInactive();
1625 template<
typename FloatTreeT>
1629 (*this)(mLeafs.getRange());
1630 mDistTree.pruneInactive();
1631 mIndexTree.pruneInactive();
1635 template<
typename FloatTreeT>
1637 FloatTreeT& distTree,
1641 : mDistTree(distTree)
1643 , mIndexTree(indexTree)
1644 , mIntersectionTree(intersectionTree)
1649 template<
typename FloatTreeT>
1652 : mDistTree(rhs.mDistTree)
1653 , mLeafs(rhs.mLeafs)
1654 , mIndexTree(rhs.mIndexTree)
1655 , mIntersectionTree(rhs.mIntersectionTree)
1660 template<
typename FloatTreeT>
1663 const tbb::blocked_range<size_t>& range)
const
1670 typename DistLeafT::ValueOnCIter iter;
1671 const FloatValueT distBG = mDistTree.background();
1672 const Int32 indexBG = mIntersectionTree.background();
1679 for (
size_t n = range.begin(); n < range.end(); ++n) {
1683 ijk = distLeaf.getOrigin();
1688 iter = distLeaf.cbeginValueOn();
1689 for (; iter; ++iter) {
1691 value = iter.getValue();
1692 if(value > 0.0)
continue;
1694 offset = iter.pos();
1695 if (maskLeaf && maskLeaf->isValueOn(offset))
continue;
1697 ijk = iter.getCoord();
1699 for (
Int32 m = 0; m < 18; ++m) {
1700 m_ijk = ijk + util::COORD_OFFSETS[m];
1708 distLeaf.setValueOff(offset, distBG);
1709 indexLeaf.setValueOff(offset, indexBG);
1722 template<
typename FloatTreeT>
1729 typedef typename FloatTreeT::template ValueConverter<Int32>::Type
IntTreeT;
1732 typedef typename FloatTreeT::template ValueConverter<bool>::Type
BoolTreeT;
1740 const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList);
1745 void operator()(
const tbb::blocked_range<size_t>&);
1759 double closestPrimDist(
const Coord&, std::vector<Int32>&,
Int32&)
const;
1763 FloatTreeT& mDistTree;
1767 const FloatValueT mExteriorBandWidth, mInteriorBandWidth, mVoxelSize;
1768 const std::vector<Vec3s>& mPointList;
1769 const std::vector<Vec4I>& mPolygonList;
1771 FloatTreeT mNewDistTree;
1777 template<
typename FloatTreeT>
1780 FloatTreeT& distTree,
1786 const std::vector<Vec3s>& pointList,
1787 const std::vector<Vec4I>& polygonList)
1789 , mDistTree(distTree)
1790 , mIndexTree(indexTree)
1791 , mMaskTree(maskTree)
1792 , mExteriorBandWidth(exteriorBandWidth)
1793 , mInteriorBandWidth(interiorBandWidth)
1794 , mVoxelSize(voxelSize)
1795 , mPointList(pointList)
1796 , mPolygonList(polygonList)
1799 , mNewMaskTree(false)
1804 template<
typename FloatTreeT>
1806 : mMaskLeafs(rhs.mMaskLeafs)
1807 , mDistTree(rhs.mDistTree)
1808 , mIndexTree(rhs.mIndexTree)
1809 , mMaskTree(rhs.mMaskTree)
1810 , mExteriorBandWidth(rhs.mExteriorBandWidth)
1811 , mInteriorBandWidth(rhs.mInteriorBandWidth)
1812 , mVoxelSize(rhs.mVoxelSize)
1813 , mPointList(rhs.mPointList)
1814 , mPolygonList(rhs.mPolygonList)
1817 , mNewMaskTree(false)
1822 template<
typename FloatTreeT>
1826 tbb::parallel_reduce(mMaskLeafs.getRange(), *
this);
1828 mDistTree.merge(mNewDistTree);
1829 mIndexTree.merge(mNewIndexTree);
1832 mMaskTree.merge(mNewMaskTree);
1836 template<
typename FloatTreeT>
1840 (*this)(mMaskLeafs.getRange());
1843 mDistTree.merge(mNewDistTree);
1844 mIndexTree.merge(mNewIndexTree);
1847 mMaskTree.merge(mNewMaskTree);
1851 template<
typename FloatTreeT>
1856 Int32 closestPrim = 0;
1870 std::vector<Int32> primitives(18);
1872 for (
size_t n = range.begin(); n < range.end(); ++n) {
1874 BoolLeafT& maskLeaf = mMaskLeafs.leaf(n);
1876 if (maskLeaf.isEmpty())
continue;
1878 ijk = maskLeaf.getOrigin();
1883 newDistAcc.
addLeaf(distLeafPt);
1888 if (!indexLeafPt) indexLeafPt = newIndexAcc.
touchLeaf(ijk);
1889 IntLeafT& indexLeaf = *indexLeafPt;
1891 bbox = maskLeaf.getNodeBoundingBox();
1894 typename BoolLeafT::ValueOnIter iter = maskLeaf.beginValueOn();
1895 for (; iter; ++iter) {
1897 ijk = iter.getCoord();
1900 distance = evalVoxelDist(ijk, distLeaf, indexLeaf, maskLeaf,
1901 primitives, closestPrim);
1903 distance = evalVoxelDist(ijk, distAcc, indexAcc, maskAcc,
1904 primitives, closestPrim);
1909 inside = distLeaf.getValue(ijk) <
FloatValueT(0.0);
1911 if (!inside && distance < mExteriorBandWidth) {
1912 distLeaf.setValueOn(pos, distance);
1913 indexLeaf.setValueOn(pos, closestPrim);
1914 }
else if (inside && distance < mInteriorBandWidth) {
1915 distLeaf.setValueOn(pos, -distance);
1916 indexLeaf.setValueOn(pos, closestPrim);
1921 for (
Int32 i = 0; i < 6; ++i) {
1922 newMaskAcc.
setValueOn(ijk + util::COORD_OFFSETS[i]);
1929 template<
typename FloatTreeT>
1933 FloatAccessorT& distAcc,
1934 IntAccessorT& indexAcc,
1935 BoolAccessorT& maskAcc,
1936 std::vector<Int32>& prims,
1937 Int32& closestPrim)
const
1944 for (
Int32 n = 0; n < 18; ++n) {
1945 n_ijk = ijk + util::COORD_OFFSETS[n];
1946 if (!maskAcc.isValueOn(n_ijk) && distAcc.probeValue(n_ijk, tmpDist)) {
1947 prims.push_back(indexAcc.getValue(n_ijk));
1948 tmpDist = std::abs(tmpDist);
1949 if (tmpDist < minDist) minDist = tmpDist;
1954 tmpDist = FloatValueT(closestPrimDist(ijk, prims, closestPrim));
1958 return tmpDist > minDist ? tmpDist : minDist + mVoxelSize;
1963 template<
typename FloatTreeT>
1965 ExpandNB<FloatTreeT>::evalVoxelDist(
1967 FloatLeafT& distLeaf,
1968 IntLeafT& indexLeaf,
1969 BoolLeafT& maskLeaf,
1970 std::vector<Int32>& prims,
1971 Int32& closestPrim)
const
1977 for (
Int32 n = 0; n < 18; ++n) {
1978 pos = FloatLeafT::coord2offset(ijk + util::COORD_OFFSETS[n]);
1979 if (!maskLeaf.isValueOn(pos) && distLeaf.probeValue(pos, tmpDist)) {
1980 prims.push_back(indexLeaf.getValue(pos));
1981 tmpDist = std::abs(tmpDist);
1982 if (tmpDist < minDist) minDist = tmpDist;
1986 tmpDist = FloatValueT(closestPrimDist(ijk, prims, closestPrim));
1987 return tmpDist > minDist ? tmpDist : minDist + mVoxelSize;
1991 template<
typename FloatTreeT>
1993 ExpandNB<FloatTreeT>::closestPrimDist(
const Coord& ijk,
1994 std::vector<Int32>& prims,
Int32& closestPrim)
const
1996 std::sort(prims.begin(), prims.end());
1998 Int32 lastPrim = -1;
1999 Vec3d uvw, voxelCenter(ijk[0], ijk[1], ijk[2]);
2002 for (
size_t n = 0, N = prims.size(); n < N; ++n) {
2003 if (prims[n] == lastPrim)
continue;
2005 lastPrim = prims[n];
2007 const Vec4I& verts = mPolygonList[lastPrim];
2010 const Vec3d a(mPointList[verts[0]]);
2011 const Vec3d b(mPointList[verts[1]]);
2012 const Vec3d c(mPointList[verts[2]]);
2014 primDist = (voxelCenter -
2019 const Vec3d d(mPointList[verts[3]]);
2021 tmpDist = (voxelCenter -
2024 if (tmpDist < primDist) primDist = tmpDist;
2027 if (primDist < dist) {
2029 closestPrim = lastPrim;
2033 return std::sqrt(dist) * double(mVoxelSize);
2037 template<
typename FloatTreeT>
2041 mNewDistTree.merge(rhs.mNewDistTree);
2042 mNewIndexTree.merge(rhs.mNewIndexTree);
2043 mNewMaskTree.merge(rhs.mNewMaskTree);
2050 template<
typename ValueType>
2054 : mVoxelSize(voxelSize)
2055 , mUnsigned(unsignedDist)
2059 template <
typename LeafNodeType>
2066 typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
2067 for (; iter; ++iter) {
2068 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
2069 val = w[!mUnsigned && int(val < ValueType(0.0))] * std::sqrt(std::abs(val));
2074 ValueType mVoxelSize;
2075 const bool mUnsigned;
2079 template<
typename ValueType>
2083 : mExBandWidth(exBandWidth)
2084 , mInBandWidth(inBandWidth)
2088 template <
typename LeafNodeType>
2091 ValueType bgValues[2];
2092 bgValues[0] = mExBandWidth;
2093 bgValues[1] = -mInBandWidth;
2095 typename LeafNodeType::ValueOffIter iter = leaf.beginValueOff();
2097 for (; iter; ++iter) {
2098 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
2099 val = bgValues[int(val < ValueType(0.0))];
2104 ValueType mExBandWidth, mInBandWidth;
2108 template<
typename ValueType>
2111 TrimOp(ValueType exBandWidth, ValueType inBandWidth)
2112 : mExBandWidth(exBandWidth)
2113 , mInBandWidth(inBandWidth)
2117 template <
typename LeafNodeType>
2120 typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
2122 for (; iter; ++iter) {
2123 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
2124 const bool inside = val < ValueType(0.0);
2126 if (inside && !(val > -mInBandWidth)) {
2127 val = -mInBandWidth;
2129 }
else if (!inside && !(val < mExBandWidth)) {
2137 ValueType mExBandWidth, mInBandWidth;
2141 template<
typename ValueType>
2148 template <
typename LeafNodeType>
2151 typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
2152 for (; iter; ++iter) {
2153 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
2163 template<
typename Gr
idType,
typename ValueType>
2167 typedef typename Scheme::template ISStencil<GridType>::StencilType
Stencil;
2174 , mVoxelSize(voxelSize)
2181 template <
typename LeafNodeType>
2184 const ValueType dt = mCFL * mVoxelSize, one(1.0), invDx = one / mVoxelSize;
2188 typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
2189 for (; iter; ++iter) {
2190 stencil.moveTo(iter);
2192 const ValueType normSqGradPhi =
2195 const ValueType phi0 = iter.getValue();
2196 const ValueType diff =
math::Sqrt(normSqGradPhi) * invDx - one;
2199 buffer.setValue(iter.pos(), phi0 - dt * S * diff);
2206 ValueType mVoxelSize, mCFL;
2210 template<
typename TreeType,
typename ValueType>
2218 template <
typename LeafNodeType>
2222 typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
2224 for (; iter; ++iter) {
2225 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
2226 val =
std::min(val, buffer.getValue(iter.pos()));
2235 template<
typename TreeType,
typename ValueType>
2243 , mBufferIndex(bufferIndex)
2247 template <
typename LeafNodeType>
2251 typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
2254 for (; iter; ++iter) {
2255 offset = iter.pos();
2256 leaf.setValueOnly(offset, buffer.getValue(offset));
2262 const size_t mBufferIndex;
2266 template<
typename TreeType>
2274 template <
typename LeafNodeType>
2277 const LeafNodeT* rhsLeaf = mTree.probeConstLeaf(leaf.getOrigin());
2278 if (rhsLeaf) leaf.topologyDifference(*rhsLeaf,
false);
2294 template<
typename FloatGr
idT,
typename InterruptT>
2296 openvdb::math::Transform::Ptr& transform,
int conversionFlags,
2297 InterruptT *interrupter,
int signSweeps)
2298 : mTransform(transform)
2299 , mConversionFlags(conversionFlags)
2300 , mSignSweeps(signSweeps)
2301 , mInterrupter(interrupter)
2304 mSignSweeps =
std::min(mSignSweeps, 1);
2308 template<
typename FloatGr
idT,
typename InterruptT>
2314 mIntersectingVoxelsGrid = BoolGridT::create(
false);
2318 template<
typename FloatGr
idT,
typename InterruptT>
2321 const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList,
2327 const FloatValueT vs = mTransform->voxelSize()[0];
2328 doConvert(pointList, polygonList, vs * exBandWidth, vs * inBandWidth);
2333 template<
typename FloatGr
idT,
typename InterruptT>
2336 const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList,
2341 const FloatValueT vs = mTransform->voxelSize()[0];
2342 doConvert(pointList, polygonList, vs * exBandWidth, 0.0,
true);
2347 template<
typename FloatGr
idT,
typename InterruptT>
2350 const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList,
2351 FloatValueT exBandWidth, FloatValueT inBandWidth,
bool unsignedDistField)
2353 mDistGrid->setTransform(mTransform);
2354 mIndexGrid->setTransform(mTransform);
2365 voxelizer(pointList, polygonList, mInterrupter);
2371 mDistGrid->tree().merge(voxelizer.sqrDistTree());
2372 mIndexGrid->tree().merge(voxelizer.primIndexTree());
2373 mIntersectingVoxelsGrid->tree().merge(voxelizer.intersectionTree());
2376 if (!unsignedDistField) {
2380 internal::ContourTracer<FloatTreeT, InterruptT> trace(
2381 mDistGrid->tree(), mIntersectingVoxelsGrid->tree(), mInterrupter);
2382 for (
int i = 0; i < mSignSweeps; ++i) {
2386 trace.runParallel();
2391 BoolTreeT signMaskTree(
false);
2394 internal::SignMask<FloatTreeT, InterruptT> signMaskOp(leafs,
2395 mDistGrid->tree(), mIntersectingVoxelsGrid->tree(), mInterrupter);
2396 signMaskOp.runParallel();
2397 signMaskTree.merge(signMaskOp.signMaskTree());
2404 if(leafs.leafCount() == 0)
break;
2406 internal::PropagateSign<FloatTreeT, InterruptT> sign(leafs,
2407 mDistGrid->tree(), mIntersectingVoxelsGrid->tree(), mInterrupter);
2409 signMaskTree.clear();
2411 signMaskTree.merge(sign.signMaskTree());
2420 tree::LeafManager<BoolTreeT> leafs(mIntersectingVoxelsGrid->tree());
2423 internal::IntersectingVoxelSign<FloatTreeT> sign(pointList, polygonList,
2424 mDistGrid->tree(), mIndexGrid->tree(), mIntersectingVoxelsGrid->tree(), leafs);
2432 internal::IntersectingVoxelCleaner<FloatTreeT> cleaner(mDistGrid->tree(),
2433 mIndexGrid->tree(), mIntersectingVoxelsGrid->tree(), leafs);
2434 cleaner.runParallel();
2441 tree::LeafManager<FloatTreeT> leafs(mDistGrid->tree());
2443 internal::ShellVoxelCleaner<FloatTreeT> cleaner(mDistGrid->tree(),
2444 leafs, mIndexGrid->tree(), mIntersectingVoxelsGrid->tree());
2446 cleaner.runParallel();
2452 inBandWidth = FloatValueT(0.0);
2455 if (mDistGrid->activeVoxelCount() == 0)
return;
2457 mIntersectingVoxelsGrid->clear();
2458 const FloatValueT voxelSize(mTransform->voxelSize()[0]);
2462 tree::LeafManager<FloatTreeT> leafs(mDistGrid->tree());
2463 leafs.foreach(internal::SqrtAndScaleOp<FloatValueT>(voxelSize, unsignedDistField));
2470 if (!unsignedDistField) {
2472 mDistGrid->tree().signedFloodFill();
2477 tree::LeafManager<FloatTreeT> leafs(mDistGrid->tree());
2479 leafs.foreach(internal::VoxelSignOp<FloatValueT>(exBandWidth, inBandWidth));
2483 FloatValueT bgValues[2];
2484 bgValues[0] = exBandWidth;
2485 bgValues[1] = -inBandWidth;
2487 typename FloatTreeT::ValueAllIter tileIt(mDistGrid->tree());
2488 tileIt.setMaxDepth(FloatTreeT::ValueAllIter::LEAF_DEPTH - 1);
2490 for ( ; tileIt; ++tileIt) {
2491 FloatValueT& val =
const_cast<FloatValueT&
>(tileIt.getValue());
2492 val = bgValues[int(val < FloatValueT(0.0))];
2498 typename FloatTreeT::Ptr newTree(
new FloatTreeT(exBandWidth));
2499 newTree->merge(mDistGrid->tree());
2500 mDistGrid->setTree(newTree);
2506 const FloatValueT minWidth = voxelSize * 2.0;
2507 if (inBandWidth > minWidth || exBandWidth > minWidth) {
2510 BoolTreeT maskTree(
false);
2511 tree::ValueAccessor<BoolTreeT> acc(maskTree);
2512 maskTree.topologyUnion(mDistGrid->tree());
2516 internal::LeafTopologyDiffOp<FloatTreeT> diffOp(mDistGrid->tree());
2521 float progress = 48, step = 0.0;
2524 2.0 * std::ceil((
std::max(inBandWidth, exBandWidth) - minWidth) / voxelSize);
2525 if (estimated <
double(maxIterations)) {
2526 maxIterations = unsigned(estimated);
2527 step = 42.0 / float(maxIterations);
2535 tree::LeafManager<BoolTreeT> leafs(maskTree);
2537 if (leafs.leafCount() == 0)
break;
2539 leafs.foreach(diffOp);
2541 internal::ExpandNB<FloatTreeT> expand(
2542 leafs, mDistGrid->tree(), mIndexGrid->tree(), maskTree,
2543 exBandWidth, inBandWidth, voxelSize, pointList, polygonList);
2545 expand.runParallel();
2547 if ((++count) >= maxIterations)
break;
2559 if (!unsignedDistField) {
2561 mDistGrid->tree().pruneLevelSet();
2562 tree::LeafManager<FloatTreeT> leafs(mDistGrid->tree(), 1);
2564 const FloatValueT offset = 0.8 * voxelSize;
2567 internal::OffsetOp<FloatValueT> offsetOp(-offset);
2569 leafs.foreach(offsetOp);
2573 leafs.foreach(internal::RenormOp<FloatGridT, FloatValueT>(*mDistGrid, leafs, voxelSize));
2575 leafs.foreach(internal::MinOp<FloatTreeT, FloatValueT>(leafs));
2579 offsetOp.resetOffset(offset - internal::Tolerance<FloatValueT>::epsilon());
2580 leafs.foreach(offsetOp);
2585 const FloatValueT minTrimWidth = voxelSize * 4.0;
2586 if (inBandWidth < minTrimWidth || exBandWidth < minTrimWidth) {
2592 tree::LeafManager<FloatTreeT> leafs(mDistGrid->tree());
2593 leafs.foreach(internal::TrimOp<FloatValueT>(
2594 exBandWidth, unsignedDistField ? exBandWidth : inBandWidth));
2600 mDistGrid->tree().pruneLevelSet();
2601 mDistGrid->tree().signedFloodFill();
2609 template<
typename Gr
idType>
2610 inline typename boost::enable_if<boost::is_floating_point<typename GridType::ValueType>,
2611 typename GridType::Ptr>::type
2613 const openvdb::math::Transform& xform,
2614 const std::vector<Vec3s>& points,
2615 const std::vector<Vec3I>& triangles,
2616 const std::vector<Vec4I>& quads,
2619 bool unsignedDistanceField =
false)
2621 std::vector<Vec3s> indexSpacePoints(points.size());
2629 std::vector<Vec4I> primitives(triangles.size() + quads.size());
2631 for (
size_t n = 0, N = triangles.size(); n < N; ++n) {
2632 Vec4I& prim = primitives[n];
2633 const Vec3I& triangle = triangles[n];
2634 prim[0] = triangle[0];
2635 prim[1] = triangle[1];
2636 prim[2] = triangle[2];
2640 for (
size_t n = 0, N = quads.size(); n < N; ++n) {
2641 primitives[n + triangles.size()] = quads[n];
2644 typename GridType::ValueType exWidth(exBandWidth);
2645 typename GridType::ValueType inWidth(inBandWidth);
2651 if (!unsignedDistanceField) {
2663 template<
typename Gr
idType>
2664 inline typename boost::disable_if<boost::is_floating_point<typename GridType::ValueType>,
2665 typename GridType::Ptr>::type
2668 const std::vector<Vec3s>& ,
2669 const std::vector<Vec3I>& ,
2670 const std::vector<Vec4I>& ,
2673 bool unsignedDistanceField =
false)
2676 "mesh to volume conversion is supported only for scalar, floating-point grids");
2683 template<
typename Gr
idType>
2684 inline typename GridType::Ptr
2686 const openvdb::math::Transform& xform,
2687 const std::vector<Vec3s>& points,
2688 const std::vector<Vec3I>& triangles,
2691 std::vector<Vec4I> quads(0);
2692 return doMeshConversion<GridType>(xform, points, triangles, quads,
2693 halfWidth, halfWidth);
2697 template<
typename Gr
idType>
2698 inline typename GridType::Ptr
2700 const openvdb::math::Transform& xform,
2701 const std::vector<Vec3s>& points,
2702 const std::vector<Vec4I>& quads,
2705 std::vector<Vec3I> triangles(0);
2706 return doMeshConversion<GridType>(xform, points, triangles, quads,
2707 halfWidth, halfWidth);
2711 template<
typename Gr
idType>
2712 inline typename GridType::Ptr
2714 const openvdb::math::Transform& xform,
2715 const std::vector<Vec3s>& points,
2716 const std::vector<Vec3I>& triangles,
2717 const std::vector<Vec4I>& quads,
2720 return doMeshConversion<GridType>(xform, points, triangles, quads,
2721 halfWidth, halfWidth);
2725 template<
typename Gr
idType>
2726 inline typename GridType::Ptr
2728 const openvdb::math::Transform& xform,
2729 const std::vector<Vec3s>& points,
2730 const std::vector<Vec3I>& triangles,
2731 const std::vector<Vec4I>& quads,
2735 return doMeshConversion<GridType>(xform, points, triangles,
2736 quads, exBandWidth, inBandWidth);
2740 template<
typename Gr
idType>
2741 inline typename GridType::Ptr
2743 const openvdb::math::Transform& xform,
2744 const std::vector<Vec3s>& points,
2745 const std::vector<Vec3I>& triangles,
2746 const std::vector<Vec4I>& quads,
2749 return doMeshConversion<GridType>(xform, points, triangles, quads,
2750 bandWidth, bandWidth,
true);
2758 inline std::ostream&
2761 ostr <<
"{[ " << rhs.
mXPrim <<
", " << rhs.
mXDist <<
"]";
2762 ostr <<
" [ " << rhs.
mYPrim <<
", " << rhs.
mYDist <<
"]";
2763 ostr <<
" [ " << rhs.
mZPrim <<
", " << rhs.
mZDist <<
"]}";
2768 inline MeshToVoxelEdgeData::EdgeData
2783 const std::vector<Vec3s>& pointList,
2784 const std::vector<Vec4I>& polygonList);
2786 void run(
bool threaded =
true);
2789 inline void operator() (
const tbb::blocked_range<size_t> &range);
2797 struct Primitive {
Vec3d a, b, c, d;
Int32 index; };
2799 template<
bool IsQuad>
2800 inline void voxelize(
const Primitive&);
2802 template<
bool IsQuad>
2803 inline bool evalPrimitive(
const Coord&,
const Primitive&);
2805 inline bool rayTriangleIntersection(
const Vec3d& origin,
const Vec3d& dir,
2812 const std::vector<Vec3s>& mPointList;
2813 const std::vector<Vec4I>& mPolygonList;
2817 IntTreeT mLastPrimTree;
2824 const std::vector<Vec3s>& pointList,
2825 const std::vector<Vec4I>& polygonList)
2828 , mPointList(pointList)
2829 , mPolygonList(polygonList)
2831 , mLastPrimAccessor(mLastPrimTree)
2840 , mPointList(rhs.mPointList)
2841 , mPolygonList(rhs.mPolygonList)
2843 , mLastPrimAccessor(mLastPrimTree)
2852 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPolygonList.size()), *
this);
2854 (*this)(tbb::blocked_range<size_t>(0, mPolygonList.size()));
2864 typedef boost::mpl::at<NodeTypeVec, boost::mpl::int_<1> >::type InternalNodeType;
2872 for ( ; leafIt; ++leafIt) {
2873 ijk = leafIt->getOrigin();
2879 mAccessor.addLeaf(rhs.mAccessor.
probeLeaf(ijk));
2880 InternalNodeType* node = rhs.mAccessor.
getNode<InternalNodeType>();
2882 rhs.mAccessor.
clear();
2892 if (!lhsLeafPt->isValueOn(offset)) {
2893 lhsLeafPt->setValueOn(offset, rhsValue);
2925 for (
size_t n = range.begin(); n < range.end(); ++n) {
2927 const Vec4I& verts = mPolygonList[n];
2929 prim.index =
Int32(n);
2930 prim.a =
Vec3d(mPointList[verts[0]]);
2931 prim.b =
Vec3d(mPointList[verts[1]]);
2932 prim.c =
Vec3d(mPointList[verts[2]]);
2935 prim.d =
Vec3d(mPointList[verts[3]]);
2936 voxelize<true>(prim);
2938 voxelize<false>(prim);
2944 template<
bool IsQuad>
2946 MeshToVoxelEdgeData::GenEdgeData::voxelize(
const Primitive& prim)
2948 std::deque<Coord> coordList;
2952 coordList.push_back(ijk);
2954 evalPrimitive<IsQuad>(ijk, prim);
2956 while (!coordList.empty()) {
2958 ijk = coordList.back();
2959 coordList.pop_back();
2961 for (
Int32 i = 0; i < 26; ++i) {
2962 nijk = ijk + util::COORD_OFFSETS[i];
2964 if (prim.index != mLastPrimAccessor.getValue(nijk)) {
2965 mLastPrimAccessor.setValue(nijk, prim.index);
2966 if(evalPrimitive<IsQuad>(nijk, prim)) coordList.push_back(nijk);
2973 template<
bool IsQuad>
2975 MeshToVoxelEdgeData::GenEdgeData::evalPrimitive(
const Coord& ijk,
const Primitive& prim)
2977 Vec3d uvw, org(ijk[0], ijk[1], ijk[2]);
2978 bool intersecting =
false;
2982 mAccessor.probeValue(ijk, edgeData);
2985 double dist = (org -
2988 if (rayTriangleIntersection(org,
Vec3d(1.0, 0.0, 0.0), prim.a, prim.c, prim.b, t)) {
2989 if (t < edgeData.mXDist) {
2990 edgeData.mXDist = t;
2991 edgeData.mXPrim = prim.index;
2992 intersecting =
true;
2996 if (rayTriangleIntersection(org,
Vec3d(0.0, 1.0, 0.0), prim.a, prim.c, prim.b, t)) {
2997 if (t < edgeData.mYDist) {
2998 edgeData.mYDist = t;
2999 edgeData.mYPrim = prim.index;
3000 intersecting =
true;
3004 if (rayTriangleIntersection(org,
Vec3d(0.0, 0.0, 1.0), prim.a, prim.c, prim.b, t)) {
3005 if (t < edgeData.mZDist) {
3006 edgeData.mZDist = t;
3007 edgeData.mZPrim = prim.index;
3008 intersecting =
true;
3014 double secondDist = (org -
3017 if (secondDist < dist) dist = secondDist;
3019 if (rayTriangleIntersection(org,
Vec3d(1.0, 0.0, 0.0), prim.a, prim.d, prim.c, t)) {
3020 if (t < edgeData.mXDist) {
3021 edgeData.mXDist = t;
3022 edgeData.mXPrim = prim.index;
3023 intersecting =
true;
3027 if (rayTriangleIntersection(org,
Vec3d(0.0, 1.0, 0.0), prim.a, prim.d, prim.c, t)) {
3028 if (t < edgeData.mYDist) {
3029 edgeData.mYDist = t;
3030 edgeData.mYPrim = prim.index;
3031 intersecting =
true;
3035 if (rayTriangleIntersection(org,
Vec3d(0.0, 0.0, 1.0), prim.a, prim.d, prim.c, t)) {
3036 if (t < edgeData.mZDist) {
3037 edgeData.mZDist = t;
3038 edgeData.mZPrim = prim.index;
3039 intersecting =
true;
3044 if (intersecting) mAccessor.setValue(ijk, edgeData);
3046 return (dist < 0.86602540378443861);
3051 MeshToVoxelEdgeData::GenEdgeData::rayTriangleIntersection(
3062 double divisor = s1.
dot(e1);
3063 if (!(std::abs(divisor) > 0.0))
return false;
3067 double inv_divisor = 1.0 / divisor;
3068 Vec3d d = origin - a;
3069 double b1 = d.
dot(s1) * inv_divisor;
3071 if (b1 < 0.0 || b1 > 1.0)
return false;
3074 double b2 = dir.
dot(s2) * inv_divisor;
3076 if (b2 < 0.0 || (b1 + b2) > 1.0)
return false;
3080 t = e2.dot(s2) * inv_divisor;
3081 return (t < 0.0) ?
false :
true;
3097 const std::vector<Vec3s>& pointList,
3098 const std::vector<Vec4I>& polygonList)
3112 std::vector<Vec3d>& points,
3113 std::vector<Index32>& primitives)
3123 point[0] = double(coord[0]) + data.
mXDist;
3124 point[1] = double(coord[1]);
3125 point[2] = double(coord[2]);
3127 points.push_back(point);
3128 primitives.push_back(data.
mXPrim);
3132 point[0] = double(coord[0]);
3133 point[1] = double(coord[1]) + data.
mYDist;
3134 point[2] = double(coord[2]);
3136 points.push_back(point);
3137 primitives.push_back(data.
mYPrim);
3141 point[0] = double(coord[0]);
3142 point[1] = double(coord[1]);
3143 point[2] = double(coord[2]) + data.
mZDist;
3145 points.push_back(point);
3146 primitives.push_back(data.
mZPrim);
3156 point[0] = double(coord[0]);
3157 point[1] = double(coord[1]) + data.
mYDist;
3158 point[2] = double(coord[2]);
3160 points.push_back(point);
3161 primitives.push_back(data.
mYPrim);
3165 point[0] = double(coord[0]);
3166 point[1] = double(coord[1]);
3167 point[2] = double(coord[2]) + data.
mZDist;
3169 points.push_back(point);
3170 primitives.push_back(data.
mZPrim);
3178 point[0] = double(coord[0]);
3179 point[1] = double(coord[1]) + data.
mYDist;
3180 point[2] = double(coord[2]);
3182 points.push_back(point);
3183 primitives.push_back(data.
mYPrim);
3192 point[0] = double(coord[0]) + data.
mXDist;
3193 point[1] = double(coord[1]);
3194 point[2] = double(coord[2]);
3196 points.push_back(point);
3197 primitives.push_back(data.
mXPrim);
3201 point[0] = double(coord[0]);
3202 point[1] = double(coord[1]) + data.
mYDist;
3203 point[2] = double(coord[2]);
3205 points.push_back(point);
3206 primitives.push_back(data.
mYPrim);
3216 point[0] = double(coord[0]) + data.
mXDist;
3217 point[1] = double(coord[1]);
3218 point[2] = double(coord[2]);
3220 points.push_back(point);
3221 primitives.push_back(data.
mXPrim);
3230 point[0] = double(coord[0]) + data.
mXDist;
3231 point[1] = double(coord[1]);
3232 point[2] = double(coord[2]);
3234 points.push_back(point);
3235 primitives.push_back(data.
mXPrim);
3239 point[0] = double(coord[0]);
3240 point[1] = double(coord[1]);
3241 point[2] = double(coord[2]) + data.
mZDist;
3243 points.push_back(point);
3244 primitives.push_back(data.
mZPrim);
3253 point[0] = double(coord[0]);
3254 point[1] = double(coord[1]);
3255 point[2] = double(coord[2]) + data.
mZDist;
3257 points.push_back(point);
3258 primitives.push_back(data.
mZPrim);
3268 #endif // OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED