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(
90 const std::vector<Vec3s>& pointList,
91 const std::vector<Vec4I>& polygonList,
103 void convertToUnsignedDistanceField(
const std::vector<Vec3s>& pointList,
104 const std::vector<Vec4I>& polygonList,
DistValueT exBandWidth);
119 void doConvert(
const std::vector<Vec3s>&,
const std::vector<Vec4I>&,
120 DistValueT exBandWidth, DistValueT inBandWidth,
bool unsignedDistField =
false);
122 openvdb::math::Transform::Ptr mTransform;
123 int mConversionFlags, mSignSweeps;
125 typename DistGridT::Ptr mDistGrid;
126 typename IndexGridT::Ptr mIndexGrid;
127 typename StencilGridT::Ptr mIntersectingVoxelsGrid;
129 InterruptT *mInterrupter;
142 template<
typename ValueType>
145 static ValueType
epsilon() {
return ValueType(1e-7); }
150 template<
typename DistTreeT,
typename IndexTreeT>
152 combine(DistTreeT& lhsDist, IndexTreeT& lhsIndex, DistTreeT& rhsDist, IndexTreeT& rhsIndex)
154 typedef typename DistTreeT::ValueType DistValueT;
158 typename DistTreeT::LeafCIter iter = rhsDist.cbeginLeaf();
163 for ( ; iter; ++iter) {
164 typename DistTreeT::LeafNodeType::ValueOnCIter it = iter->cbeginValueOn();
169 rhsValue = it.getValue();
170 DistValueT& lhsValue =
const_cast<DistValueT&
>(lhsDistAccessor.
getValue(ijk));
172 if (-rhsValue < std::abs(lhsValue)) {
191 template<
typename DistTreeT,
typename InterruptT = util::NullInterrupter>
199 typedef typename DistTreeT::template ValueConverter<Int32>::Type
IndexTreeT;
201 typedef typename DistTreeT::template ValueConverter<bool>::Type
StencilTreeT;
206 const std::vector<Vec4I>& polygonList, InterruptT *interrupter = NULL);
214 void operator() (
const tbb::blocked_range<size_t> &range);
226 bool evalVoxel(
const Coord& ijk,
const Int32 polyIdx);
228 std::vector<Vec3s>
const *
const mPointList;
229 std::vector<Vec4I>
const *
const mPolygonList;
231 DistTreeT mSqrDistTree;
232 DistAccessorT mSqrDistAccessor;
234 IndexTreeT mPrimIndexTree;
235 IndexAccessorT mPrimIndexAccessor;
237 StencilTreeT mIntersectionTree;
238 StencilAccessorT mIntersectionAccessor;
241 IndexTreeT mLastPrimTree;
242 IndexAccessorT mLastPrimAccessor;
244 InterruptT *mInterrupter;
248 template<
typename DistTreeT,
typename InterruptT>
252 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPolygonList->size()), *
this);
255 template<
typename DistTreeT,
typename InterruptT>
259 (*this)(tbb::blocked_range<size_t>(0, mPolygonList->size()));
262 template<
typename DistTreeT,
typename InterruptT>
264 const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList,
265 InterruptT *interrupter)
266 : mPointList(&pointList)
267 , mPolygonList(&polygonList)
269 , mSqrDistAccessor(mSqrDistTree)
271 , mPrimIndexAccessor(mPrimIndexTree)
272 , mIntersectionTree(false)
273 , mIntersectionAccessor(mIntersectionTree)
275 , mLastPrimAccessor(mLastPrimTree)
276 , mInterrupter(interrupter)
280 template<
typename DistTreeT,
typename InterruptT>
283 : mPointList(rhs.mPointList)
284 , mPolygonList(rhs.mPolygonList)
286 , mSqrDistAccessor(mSqrDistTree)
288 , mPrimIndexAccessor(mPrimIndexTree)
289 , mIntersectionTree(false)
290 , mIntersectionAccessor(mIntersectionTree)
292 , mLastPrimAccessor(mLastPrimTree)
293 , mInterrupter(rhs.mInterrupter)
297 template<
typename DistTreeT,
typename InterruptT>
302 double edge_max = std::abs(v1[0] - v0[0]);
303 edge_max =
std::max(edge_max, std::abs(v1[1] - v0[1]));
304 edge_max =
std::max(edge_max, std::abs(v1[2] - v0[2]));
305 edge_max =
std::max(edge_max, std::abs(v0[0] - v2[0]));
306 edge_max =
std::max(edge_max, std::abs(v0[1] - v2[1]));
307 edge_max =
std::max(edge_max, std::abs(v0[2] - v2[2]));
308 return edge_max < 200.0;
311 template<
typename DistTreeT,
typename InterruptT>
315 std::deque<Coord> coordList;
320 for (
size_t n = range.begin(); n < range.end(); ++n) {
322 if (mInterrupter && mInterrupter->wasInterrupted()) {
323 tbb::task::self().cancel_group_execution();
327 const Int32 primIdx = n;
328 const Vec4I verts = (*mPolygonList)[n];
330 Vec3d p0((*mPointList)[verts[0]]);
331 Vec3d p1((*mPointList)[verts[1]]);
332 Vec3d p2((*mPointList)[verts[2]]);
334 if (shortEdge(p0, p1, p2)) {
339 evalVoxel(ijk, primIdx);
340 coordList.push_back(ijk);
344 evalVoxel(ijk, primIdx);
345 coordList.push_back(ijk);
349 evalVoxel(ijk, primIdx);
350 coordList.push_back(ijk);
353 Vec3d p3((*mPointList)[verts[3]]);
355 evalVoxel(ijk, primIdx);
356 coordList.push_back(ijk);
359 while (!coordList.empty()) {
360 if (mInterrupter && mInterrupter->wasInterrupted()) {
364 ijk = coordList.back();
365 coordList.pop_back();
367 mIntersectionAccessor.setActiveState(ijk,
true);
369 for (
Int32 i = 0; i < 26; ++i) {
372 if (primIdx != mLastPrimAccessor.getValue(n_ijk)) {
373 mLastPrimAccessor.setValue(n_ijk, n);
374 if(evalVoxel(n_ijk, n)) coordList.push_back(n_ijk);
382 evalVoxel(ijk, primIdx);
384 mLastPrimAccessor.setValue(ijk, primIdx);
387 while (!auxTree.empty()) {
389 if (mInterrupter && mInterrupter->wasInterrupted()) {
393 typename StencilTreeT::LeafIter leafIter = auxTree.beginLeaf();
394 for (; leafIter; leafIter.next()) {
396 if (mInterrupter && mInterrupter->wasInterrupted()) {
400 typename StencilTreeT::LeafNodeType::ValueOnIter iter = leafIter->beginValueOn();
401 for (; iter; iter.next()) {
402 ijk = iter.getCoord();
405 mIntersectionAccessor.setActiveState(ijk,
true);
407 for (
Int32 i = 0; i < 26; ++i) {
410 if (primIdx != mLastPrimAccessor.getValue(n_ijk)) {
411 mLastPrimAccessor.setValue(n_ijk, n);
418 auxTree.pruneInactive();
424 template<
typename DistTreeT,
typename InterruptT>
428 Vec3d voxelCenter(ijk[0], ijk[1], ijk[2]);
429 const Vec4I& verts = (*mPolygonList)[polyIdx];
430 const std::vector<Vec3s>& points = *mPointList;
433 const Vec3d a(points[verts[0]]);
434 const Vec3d b(points[verts[1]]);
435 const Vec3d c(points[verts[2]]);
438 double dist = (voxelCenter -
443 const Vec3d d(points[verts[3]]);
445 double secondDist = (voxelCenter -
448 if (secondDist < dist) dist = secondDist;
451 const DistValueT tmp(dist);
452 if (tmp < std::abs(mSqrDistAccessor.getValue(ijk))) {
453 mSqrDistAccessor.setValue(ijk, -tmp);
454 mPrimIndexAccessor.setValue(ijk, polyIdx);
457 return (dist < 0.86602540378443861);
460 template<
typename DistTreeT,
typename InterruptT>
464 typename DistTreeT::LeafCIter iter = rhs.mSqrDistTree.cbeginLeaf();
468 for ( ; iter; ++iter) {
469 typename DistTreeT::LeafNodeType::ValueOnCIter it = iter->cbeginValueOn();
474 rhsDist = it.getValue();
475 DistValueT lhsDist = mSqrDistAccessor.getValue(ijk);
477 if (-rhsDist < std::abs(lhsDist)) {
478 mSqrDistAccessor.setValue(ijk, rhsDist);
479 mPrimIndexAccessor.setValue(ijk, rhs.mPrimIndexAccessor.getValue(ijk));
484 mIntersectionTree.merge(rhs.mIntersectionTree);
494 template<
typename DistTreeT,
typename InterruptT = util::NullInterrupter>
502 typedef typename DistTreeT::template ValueConverter<bool>::Type
StencilTreeT;
513 void operator()(
const tbb::blocked_range<int> &range)
const;
518 int sparseScan(
int slice)
const;
520 DistTreeT& mDistTree;
529 std::vector<Index> mStepSize;
531 InterruptT *mInterrupter;
534 template<
typename DistTreeT,
typename InterruptT>
538 tbb::parallel_for(tbb::blocked_range<int>(mBBox.min()[0], mBBox.max()[0]+1), *
this);
541 template<
typename DistTreeT,
typename InterruptT>
545 (*this)(tbb::blocked_range<int>(mBBox.min()[0], mBBox.max()[0]+1));
548 template<
typename DistTreeT,
typename InterruptT>
550 DistTreeT& distTree,
const StencilTreeT& intersectionTree, InterruptT *interrupter)
551 : mDistTree(distTree)
552 , mDistAccessor(mDistTree)
553 , mIntersectionTree(intersectionTree)
554 , mIntersectionAccessor(mIntersectionTree)
557 , mInterrupter(interrupter)
560 std::vector<Index> dims;
561 mDistTree.getNodeLog2Dims(dims);
563 mStepSize.resize(dims.size()+1, 1);
565 for (
int idx = static_cast<int>(dims.size()) - 1; idx > -1; --idx) {
566 exponent += dims[idx];
567 mStepSize[idx] = 1 << exponent;
570 mDistTree.evalLeafBoundingBox(mBBox);
573 const int tileDim = mStepSize[0];
575 for (
size_t i = 0; i < 3; ++i) {
578 double diff = std::abs(
double(mBBox.
min()[i])) / double(tileDim);
580 if (mBBox.
min()[i] <= tileDim) {
581 n = int(std::ceil(diff));
582 mBBox.
min()[i] = - n * tileDim;
584 n = int(std::floor(diff));
585 mBBox.
min()[i] = n * tileDim;
588 n = int(std::ceil(std::abs(
double(mBBox.
max()[i] - mBBox.
min()[i])) / double(tileDim)));
589 mBBox.
max()[i] = mBBox.
min()[i] + n * tileDim;
593 template<
typename DistTreeT,
typename InterruptT>
596 : mDistTree(rhs.mDistTree)
597 , mDistAccessor(mDistTree)
598 , mIntersectionTree(rhs.mIntersectionTree)
599 , mIntersectionAccessor(mIntersectionTree)
601 , mStepSize(rhs.mStepSize)
602 , mInterrupter(rhs.mInterrupter)
606 template<
typename DistTreeT,
typename InterruptT>
612 for (
int n = range.begin(); n < range.end(); n += iStep) {
614 if (mInterrupter && mInterrupter->wasInterrupted()) {
615 tbb::task::self().cancel_group_execution();
619 iStep = sparseScan(n);
623 template<
typename DistTreeT,
typename InterruptT>
627 bool lastVoxelWasOut =
true;
630 Coord ijk(slice, mBBox.min()[1], mBBox.min()[2]);
631 Coord step(mStepSize[mDistAccessor.getValueDepth(ijk) + 1]);
634 for (ijk[1] = mBBox.
min()[1]; ijk[1] <= mBBox.
max()[1]; ijk[1] += step[1]) {
636 if (mInterrupter && mInterrupter->wasInterrupted()) {
640 step[1] = mStepSize[mDistAccessor.getValueDepth(ijk) + 1];
641 step[0] =
std::min(step[0], step[1]);
643 for (ijk[2] = mBBox.
min()[2]; ijk[2] <= mBBox.
max()[2]; ijk[2] += step[2]) {
645 step[2] = mStepSize[mDistAccessor.getValueDepth(ijk) + 1];
646 step[1] =
std::min(step[1], step[2]);
647 step[0] =
std::min(step[0], step[2]);
650 if (mDistAccessor.isValueOn(ijk)) {
653 if (mIntersectionAccessor.isValueOn(ijk)) {
655 lastVoxelWasOut =
false;
658 }
else if (lastVoxelWasOut) {
660 DistValueT& val =
const_cast<DistValueT&
>(mDistAccessor.getValue(ijk));
666 for (
Int32 n = 3; n < 6; n += 2) {
669 if (mDistAccessor.probeValue(n_ijk, val) && val > 0) {
670 lastVoxelWasOut =
true;
675 if (lastVoxelWasOut) {
677 DistValueT& v =
const_cast<DistValueT&
>(mDistAccessor.getValue(ijk));
680 const int tmp_k = ijk[2];
683 for (--ijk[2]; ijk[2] >= last_k; --ijk[2]) {
684 if (mIntersectionAccessor.isValueOn(ijk))
break;
685 DistValueT& vb =
const_cast<DistValueT&
>(mDistAccessor.getValue(ijk));
686 if(vb < DistValueT(0.0)) vb = -vb;
712 template<
typename DistTreeT>
720 typedef typename DistTreeT::template ValueConverter<Int32>::Type
IndexTreeT;
722 typedef typename DistTreeT::template ValueConverter<bool>::Type
StencilTreeT;
728 const std::vector<Vec3s>& pointList,
729 const std::vector<Vec4I>& polygonList,
741 void operator()(
const tbb::blocked_range<size_t>&)
const;
746 void evalVoxel(
const Coord& ijk)
const;
747 Vec3d getClosestPointDir(
const Coord& ijk)
const;
749 std::vector<Vec3s>
const *
const mPointList;
750 std::vector<Vec4I>
const *
const mPolygonList;
752 DistTreeT& mDistTree;
763 template<
typename DistTreeT>
767 tbb::parallel_for(mLeafs.getRange(), *
this);
770 template<
typename DistTreeT>
774 (*this)(mLeafs.getRange());
777 template<
typename DistTreeT>
779 const std::vector<Vec3s>& pointList,
780 const std::vector<Vec4I>& polygonList,
785 : mPointList(&pointList)
786 , mPolygonList(&polygonList)
787 , mDistTree(distTree)
788 , mDistAccessor(mDistTree)
789 , mIndexTree(indexTree)
790 , mIndexAccessor(mIndexTree)
791 , mIntersectionTree(intersectionTree)
792 , mIntersectionAccessor(mIntersectionTree)
797 template<
typename DistTreeT>
800 : mPointList(rhs.mPointList)
801 , mPolygonList(rhs.mPolygonList)
802 , mDistTree(rhs.mDistTree)
803 , mDistAccessor(mDistTree)
804 , mIndexTree(rhs.mIndexTree)
805 , mIndexAccessor(mIndexTree)
806 , mIntersectionTree(rhs.mIntersectionTree)
807 , mIntersectionAccessor(mIntersectionTree)
812 template<
typename DistTreeT>
815 const tbb::blocked_range<size_t>& range)
const
817 typename StencilTreeT::LeafNodeType::ValueOnCIter iter;
819 for (
size_t n = range.begin(); n < range.end(); ++n) {
820 iter = mLeafs.leaf(n).cbeginValueOn();
821 for (; iter; ++iter) {
822 evalVoxel(iter.getCoord());
827 template<
typename DistTreeT>
831 const DistValueT val = mDistAccessor.getValue(ijk),
zeroVal(0.0);
835 Vec3d dir = getClosestPointDir(ijk), n_dir;
840 for (
Int32 n = 0; n < 26; ++n) {
841 n_ijk = ijk + util::COORD_OFFSETS[n];
843 if (mIntersectionAccessor.isValueOn(n_ijk))
continue;
844 if (!mDistAccessor.probeValue(n_ijk, n_val))
continue;
847 n_dir = getClosestPointDir(n_ijk);
849 if (n_dir.dot(dir) > 0.0 ) {
851 mDistAccessor.setValue(ijk, -val);
857 template<
typename DistTreeT>
859 IntersectingVoxelSign<DistTreeT>::getClosestPointDir(
const Coord& ijk)
const
861 Vec3d voxelCenter(ijk[0], ijk[1], ijk[2]);
862 const Vec4I& prim = (*mPolygonList)[mIndexAccessor.getValue(ijk)];
863 const std::vector<Vec3s>& points = *mPointList;
866 const Vec3d a(points[prim[0]]);
867 const Vec3d b(points[prim[1]]);
868 const Vec3d c(points[prim[2]]);
876 const Vec3d d(points[prim[3]]);
897 template<
typename DistTreeT>
905 typedef typename DistTreeT::template ValueConverter<Int32>::Type
IndexTreeT;
907 typedef typename DistTreeT::template ValueConverter<bool>::Type
StencilTreeT;
921 void operator()(
const tbb::blocked_range<size_t>&)
const;
926 DistTreeT& mDistTree;
937 template<
typename DistTreeT>
941 tbb::parallel_for(mLeafs.getRange(), *
this);
942 mIntersectionTree.pruneInactive();
945 template<
typename DistTreeT>
949 (*this)(mLeafs.getRange());
950 mIntersectionTree.pruneInactive();
953 template<
typename DistTreeT>
959 : mDistTree(distTree)
960 , mDistAccessor(mDistTree)
961 , mIndexTree(indexTree)
962 , mIndexAccessor(mIndexTree)
963 , mIntersectionTree(intersectionTree)
964 , mIntersectionAccessor(mIntersectionTree)
969 template<
typename DistTreeT>
972 : mDistTree(rhs.mDistTree)
973 , mDistAccessor(mDistTree)
974 , mIndexTree(rhs.mIndexTree)
975 , mIndexAccessor(mIndexTree)
976 , mIntersectionTree(rhs.mIntersectionTree)
977 , mIntersectionAccessor(mIntersectionTree)
982 template<
typename DistTreeT>
985 const tbb::blocked_range<size_t>& range)
const
989 DistValueT value, bg = mDistTree.background();
991 typename StencilTreeT::LeafNodeType::ValueOnCIter iter;
993 for (
size_t n = range.begin(); n < range.end(); ++n) {
995 typename StencilTreeT::LeafNodeType& leaf = mLeafs.leaf(n);
996 iter = leaf.cbeginValueOn();
998 for (; iter; ++iter) {
1000 ijk = iter.getCoord();
1003 for (
Int32 m = 0; m < 26; ++m) {
1004 m_ijk = ijk + util::COORD_OFFSETS[m];
1005 if (mDistAccessor.probeValue(m_ijk, value)) {
1013 if (turnOff) leaf.setValueOff(ijk, bg);
1025 template<
typename DistTreeT>
1034 typedef typename DistTreeT::template ValueConverter<Int32>::Type
IndexTreeT;
1036 typedef typename DistTreeT::template ValueConverter<bool>::Type
StencilTreeT;
1049 void operator()(
const tbb::blocked_range<size_t>&)
const;
1054 DistTreeT& mDistTree;
1065 template<
typename DistTreeT>
1069 tbb::parallel_for(mLeafs.getRange(), *
this);
1070 mDistTree.pruneInactive();
1071 mIndexTree.pruneInactive();
1074 template<
typename DistTreeT>
1078 (*this)(mLeafs.getRange());
1079 mDistTree.pruneInactive();
1080 mIndexTree.pruneInactive();
1083 template<
typename DistTreeT>
1085 DistTreeT& distTree,
1089 : mDistTree(distTree)
1091 , mDistAccessor(mDistTree)
1092 , mIndexTree(indexTree)
1093 , mIndexAccessor(mIndexTree)
1094 , mIntersectionTree(intersectionTree)
1095 , mIntersectionAccessor(mIntersectionTree)
1099 template<
typename DistTreeT>
1102 : mDistTree(rhs.mDistTree)
1103 , mLeafs(rhs.mLeafs)
1104 , mDistAccessor(mDistTree)
1105 , mIndexTree(rhs.mIndexTree)
1106 , mIndexAccessor(mIndexTree)
1107 , mIntersectionTree(rhs.mIntersectionTree)
1108 , mIntersectionAccessor(mIntersectionTree)
1112 template<
typename DistTreeT>
1115 const tbb::blocked_range<size_t>& range)
const
1122 const DistValueT distC = -0.86602540378443861;
1123 const DistValueT distBG = mDistTree.background();
1124 const Int32 indexBG = mIntersectionTree.background();
1126 typename DistTreeT::LeafNodeType::ValueOnCIter iter;
1127 for (
size_t n = range.begin(); n < range.end(); ++n) {
1129 typename DistTreeT::LeafNodeType& leaf = mLeafs.leaf(n);
1130 iter = leaf.cbeginValueOn();
1131 for (; iter; ++iter) {
1133 value = iter.getValue();
1134 if(value > 0.0)
continue;
1136 ijk = iter.getCoord();
1137 if (mIntersectionAccessor.isValueOn(ijk))
continue;
1140 for (
Int32 m = 0; m < 18; ++m) {
1141 m_ijk = ijk + util::COORD_OFFSETS[m];
1143 if (mIntersectionAccessor.isValueOn(m_ijk)) {
1150 leaf.setValueOff(ijk, distBG);
1153 mIndexAccessor.setValueOff(ijk, indexBG);
1156 if (value > distC) leaf.setValue(ijk, distC);
1169 template<
typename DistTreeT>
1177 typedef typename DistTreeT::template ValueConverter<Int32>::Type
IndexTreeT;
1179 typedef typename DistTreeT::template ValueConverter<bool>::Type
StencilTreeT;
1184 ExpandNB(
const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList,
1195 void operator()(
const tbb::blocked_range<size_t>&)
const;
1200 double getDistToPrim(
const Coord& ijk,
const Int32 polyIdx)
const;
1202 std::vector<Vec3s>
const *
const mPointList;
1203 std::vector<Vec4I>
const *
const mPolygonList;
1205 DistTreeT& mDistTree;
1210 const DistValueT mExteriorBandWidth, mInteriorBandWidth, mVoxelSize;
1213 template<
typename DistTreeT>
1217 tbb::parallel_for(mLeafs.getRange(), *
this);
1218 mMaskTree.pruneInactive();
1221 template<
typename DistTreeT>
1225 (*this)(mLeafs.getRange());
1226 mMaskTree.pruneInactive();
1229 template<
typename DistTreeT>
1231 const std::vector<Vec3s>& pointList,
1232 const std::vector<Vec4I>& polygonList,
1233 DistTreeT& distTree,
1239 : mPointList(&pointList)
1240 , mPolygonList(&polygonList)
1241 , mDistTree(distTree)
1242 , mIndexTree(indexTree)
1243 , mMaskTree(maskTree)
1245 , mExteriorBandWidth(exteriorBandWidth)
1246 , mInteriorBandWidth(interiorBandWidth)
1247 , mVoxelSize(voxelSize)
1251 template<
typename DistTreeT>
1253 : mPointList(rhs.mPointList)
1254 , mPolygonList(rhs.mPolygonList)
1255 , mDistTree(rhs.mDistTree)
1256 , mIndexTree(rhs.mIndexTree)
1257 , mMaskTree(rhs.mMaskTree)
1258 , mLeafs(rhs.mLeafs)
1259 , mExteriorBandWidth(rhs.mExteriorBandWidth)
1260 , mInteriorBandWidth(rhs.mInteriorBandWidth)
1261 , mVoxelSize(rhs.mVoxelSize)
1265 template<
typename DistTreeT>
1269 typedef typename DistTreeT::LeafNodeType DistLeafT;
1270 typedef typename IndexTreeT::LeafNodeType IndexLeafT;
1271 typedef typename StencilTreeT::LeafNodeType StencilLeafT;
1274 Int32 closestPrimIndex = 0;
1279 for (
size_t n = range.begin(); n < range.end(); ++n) {
1281 StencilLeafT& maskLeaf = mLeafs.leaf(n);
1283 const Coord origin = maskLeaf.getOrigin();
1285 DistLeafT* distLeafPt = distAcc.
probeLeaf(origin);
1286 IndexLeafT* indexLeafPt = indexAcc.
probeLeaf(origin);
1288 if (distLeafPt != NULL && indexLeafPt != NULL) {
1290 DistLeafT& distLeaf = *distLeafPt;
1291 IndexLeafT& indexLeaf = *indexLeafPt;
1294 typename StencilLeafT::ValueOnIter iter = maskLeaf.beginValueOn();
1295 for (; iter; ++iter) {
1297 const Index pos = iter.pos();
1299 if (distLeaf.isValueOn(pos)) {
1305 getDist(iter.getCoord(), distAcc, indexAcc, closestPrimIndex);
1307 const bool inside = distLeaf.getValue(pos) <
DistValueT(0.0);
1309 if (!inside && distance < mExteriorBandWidth) {
1310 distLeaf.setValueOn(pos, distance);
1311 indexLeaf.setValueOn(pos, closestPrimIndex);
1312 }
else if (inside && distance < mInteriorBandWidth) {
1313 distLeaf.setValueOn(pos, -distance);
1314 indexLeaf.setValueOn(pos, closestPrimIndex);
1321 maskLeaf.setValuesOff();
1326 template<
typename DistTreeT>
1329 IndexAccessorT& indexAcc,
Int32& primIndex)
const
1331 Vec3d voxelCenter(ijk[0], ijk[1], ijk[2]);
1336 for (
Int32 n = 0; n < 18; ++n) {
1337 n_ijk = ijk + util::COORD_OFFSETS[n];
1338 if (distAcc.probeValue(n_ijk, nDist)) {
1339 nDist = std::abs(nDist);
1342 primIndex = indexAcc.getValue(n_ijk);
1348 return getDistToPrim(ijk, primIndex);
1352 template<
typename DistTreeT>
1354 ExpandNB<DistTreeT>::getDistToPrim(
const Coord& ijk,
const Int32 polyIdx)
const
1356 Vec3d voxelCenter(ijk[0], ijk[1], ijk[2]);
1357 const Vec4I& verts = (*mPolygonList)[polyIdx];
1358 const std::vector<Vec3s>& points = *mPointList;
1361 const Vec3d a(points[verts[0]]);
1362 const Vec3d b(points[verts[1]]);
1363 const Vec3d c(points[verts[2]]);
1366 double dist = (voxelCenter -
1371 const Vec3d d(points[verts[3]]);
1373 double secondDist = (voxelCenter -
1376 if (secondDist < dist) dist = secondDist;
1379 return std::sqrt(dist) * double(mVoxelSize);
1394 template<
typename DistTreeT>
1397 typename DistTreeT::template ValueConverter<bool>::Type& intersectionTree)
1399 typedef typename DistTreeT::template ValueConverter<bool>::Type StencilTreeT;
1402 typedef typename DistTreeT::ValueType DistValueT;
1404 StencilAccessorT intrAccessor(intersectionTree);
1405 DistAccessorT distAccessor(distTree);
1407 std::deque<Coord> coordList;
1408 coordList.push_back(seed);
1411 while (!coordList.empty()) {
1412 ijk = coordList.back();
1413 coordList.pop_back();
1415 if (!distAccessor.isValueOn(ijk))
continue;
1417 DistValueT& dist =
const_cast<DistValueT&
>(distAccessor.getValue(ijk));
1418 if (!(dist < 0.0))
continue;
1421 for (
int n = 0; n < 6; ++n) {
1422 n_ijk = ijk + util::COORD_OFFSETS[n];
1424 if (!intrAccessor.isValueOn(n_ijk)) {
1425 if (distAccessor.isValueOn(n_ijk)) {
1426 if (distAccessor.getValue(n_ijk) < 0.0) {
1427 coordList.push_back(n_ijk);
1444 template<
typename DistTreeT,
typename InterruptT>
1447 typename DistTreeT::template ValueConverter<bool>::Type& intersectionTree,
1448 InterruptT *interrupter = NULL)
1450 typedef typename DistTreeT::template ValueConverter<bool>::Type StencilTreeT;
1453 typedef typename DistTreeT::ValueType DistValueT;
1455 StencilAccessorT intrAccessor(intersectionTree);
1456 DistAccessorT distAccessor(distTree);
1459 typename DistTreeT::LeafIter leafIter = distTree.beginLeaf();
1460 for (; leafIter; leafIter.next()) {
1462 if (interrupter && interrupter->wasInterrupted())
break;
1464 typename DistTreeT::LeafNodeType::ValueOnIter iter = leafIter->beginValueOn();
1465 for (; iter; iter.next()) {
1467 ijk = iter.getCoord();
1470 if (intrAccessor.isValueOn(ijk))
continue;
1472 if (iter.getValue() < 0.0) {
1473 for (
Int32 n = 0; n < 6; ++n) {
1474 n_ijk = ijk + util::COORD_OFFSETS[n];
1476 if (distAccessor.isValueOn(n_ijk) && distAccessor.getValue(n_ijk) > 0.0) {
1488 template<
typename ValueType>
1492 : mVoxelSize(voxelSize)
1493 , mUnsigned(unsignedDist)
1497 template <
typename LeafNodeType>
1504 typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
1505 for (; iter; ++iter) {
1506 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
1507 val = w[!mUnsigned && int(val < ValueType(0.0))] * std::sqrt(std::abs(val));
1512 ValueType mVoxelSize;
1513 const bool mUnsigned;
1517 template<
typename ValueType>
1521 : mExBandWidth(exBandWidth)
1522 , mInBandWidth(inBandWidth)
1526 template <
typename LeafNodeType>
1529 ValueType bgValues[2];
1530 bgValues[0] = mExBandWidth;
1531 bgValues[1] = -mInBandWidth;
1533 typename LeafNodeType::ValueOffIter iter = leaf.beginValueOff();
1535 for (; iter; ++iter) {
1536 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
1537 val = bgValues[int(val < ValueType(0.0))];
1542 ValueType mExBandWidth, mInBandWidth;
1546 template<
typename ValueType>
1549 TrimOp(ValueType exBandWidth, ValueType inBandWidth)
1550 : mExBandWidth(exBandWidth)
1551 , mInBandWidth(inBandWidth)
1555 template <
typename LeafNodeType>
1558 typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
1560 for (; iter; ++iter) {
1561 const ValueType& val = iter.getValue();
1562 const bool inside = val < ValueType(0.0);
1564 if (inside && !(val > -mInBandWidth)) {
1565 iter.setValue(-mInBandWidth);
1567 }
else if (!inside && !(val < mInBandWidth)) {
1568 iter.setValue(mExBandWidth);
1575 ValueType mExBandWidth, mInBandWidth;
1579 template<
typename ValueType>
1586 template <
typename LeafNodeType>
1589 typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
1590 for (; iter; ++iter) {
1591 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
1601 template<
typename Gr
idType,
typename ValueType>
1605 typedef typename Scheme::template ISStencil<GridType>::StencilType
Stencil;
1612 , mVoxelSize(voxelSize)
1619 template <
typename LeafNodeType>
1622 const ValueType dt = mCFL * mVoxelSize, one(1.0), invDx = one / mVoxelSize;
1627 typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
1628 for (; iter; ++iter) {
1629 stencil.moveTo(iter);
1631 const ValueType normSqGradPhi =
1634 const ValueType phi0 = stencil.getValue();
1635 const ValueType diff =
math::Sqrt(normSqGradPhi) * invDx - one;
1638 buffer.setValue(iter.pos(), phi0 - dt * S * diff);
1645 ValueType mVoxelSize, mCFL;
1649 template<
typename TreeType,
typename ValueType>
1657 template <
typename LeafNodeType>
1662 typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
1663 for (; iter; ++iter) {
1664 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
1665 val =
std::min(val, buffer.getValue(iter.pos()));
1674 template<
typename TreeType,
typename ValueType>
1682 , mBufferIndex(bufferIndex)
1686 template <
typename LeafNodeType>
1691 typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
1692 for (; iter; ++iter) {
1693 leaf.setValueOnly(iter.pos(), buffer.getValue(iter.pos()));
1699 const size_t mBufferIndex;
1710 template<
typename DistGr
idT,
typename InterruptT>
1712 openvdb::math::Transform::Ptr& transform,
int conversionFlags,
1713 InterruptT *interrupter,
int signSweeps)
1714 : mTransform(transform)
1715 , mConversionFlags(conversionFlags)
1716 , mSignSweeps(signSweeps)
1717 , mInterrupter(interrupter)
1720 mSignSweeps =
std::min(mSignSweeps, 1);
1724 template<
typename DistGr
idT,
typename InterruptT>
1730 mIntersectingVoxelsGrid = StencilGridT::create(
false);
1734 template<
typename DistGr
idT,
typename InterruptT>
1737 const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList,
1743 const DistValueT vs = mTransform->voxelSize()[0];
1744 doConvert(pointList, polygonList, vs * exBandWidth, vs * inBandWidth);
1749 template<
typename DistGr
idT,
typename InterruptT>
1752 const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList,
1757 const DistValueT vs = mTransform->voxelSize()[0];
1758 doConvert(pointList, polygonList, vs * exBandWidth, 0.0,
true);
1762 template<
typename DistGr
idT,
typename InterruptT>
1765 const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList,
1766 DistValueT exBandWidth, DistValueT inBandWidth,
bool unsignedDistField)
1768 mDistGrid->setTransform(mTransform);
1769 mIndexGrid->setTransform(mTransform);
1771 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1776 voxelizer(pointList, polygonList, mInterrupter);
1778 voxelizer.runParallel();
1780 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1782 mDistGrid->tree().merge(voxelizer.sqrDistTree());
1783 mIndexGrid->tree().merge(voxelizer.primIndexTree());
1784 mIntersectingVoxelsGrid->tree().merge(voxelizer.intersectionTree());
1787 if (!unsignedDistField) {
1792 internal::ContourTracer<DistTreeT, InterruptT> trace(
1793 mDistGrid->tree(), mIntersectingVoxelsGrid->tree(), mInterrupter);
1795 for (
int i = 0; i < mSignSweeps; ++i) {
1797 if (mInterrupter && mInterrupter->wasInterrupted())
break;
1798 trace.runParallel();
1800 if (mInterrupter && mInterrupter->wasInterrupted())
break;
1803 internal::propagateSign<DistTreeT, InterruptT>
1804 (mDistGrid->tree(), mIntersectingVoxelsGrid->tree(), mInterrupter);
1808 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1811 tree::LeafManager<StencilTreeT> leafs(mIntersectingVoxelsGrid->tree());
1814 internal::IntersectingVoxelSign<DistTreeT> sign(pointList, polygonList,
1815 mDistGrid->tree(), mIndexGrid->tree(), mIntersectingVoxelsGrid->tree(), leafs);
1819 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1823 internal::IntersectingVoxelCleaner<DistTreeT> cleaner(mDistGrid->tree(),
1824 mIndexGrid->tree(), mIntersectingVoxelsGrid->tree(), leafs);
1826 cleaner.runParallel();
1829 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1835 tree::LeafManager<DistTreeT> leafs(mDistGrid->tree());
1837 internal::ShellVoxelCleaner<DistTreeT> cleaner(mDistGrid->tree(),
1838 leafs, mIndexGrid->tree(), mIntersectingVoxelsGrid->tree());
1840 cleaner.runParallel();
1843 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1846 inBandWidth = DistValueT(0.0);
1849 if (mDistGrid->activeVoxelCount() == 0)
return;
1851 const DistValueT voxelSize(mTransform->voxelSize()[0]);
1855 typedef internal::SqrtAndScaleOp<DistValueT> XForm;
1856 tree::LeafManager<DistTreeT> leafs(mDistGrid->tree());
1857 XForm op(voxelSize, unsignedDistField);
1858 LeafTransformer<DistTreeT, XForm> transform(leafs, op);
1860 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1862 transform.runParallel();
1865 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1867 if (!unsignedDistField) {
1869 mDistGrid->tree().signedFloodFill();
1871 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1875 tree::LeafManager<DistTreeT> leafs(mDistGrid->tree());
1877 typedef internal::VoxelSignOp<DistValueT> SignXForm;
1878 SignXForm op(exBandWidth, inBandWidth);
1880 LeafTransformer<DistTreeT, SignXForm> transform(leafs, op);
1881 transform.runParallel();
1883 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1885 DistValueT bgValues[2];
1886 bgValues[0] = exBandWidth;
1887 bgValues[1] = -inBandWidth;
1889 typename DistTreeT::ValueAllIter tileIt(mDistGrid->tree());
1890 tileIt.setMaxDepth(DistTreeT::ValueAllIter::LEAF_DEPTH - 1);
1892 for ( ; tileIt; ++tileIt) {
1893 DistValueT& val =
const_cast<DistValueT&
>(tileIt.getValue());
1894 val = bgValues[int(val < DistValueT(0.0))];
1897 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1900 typename DistTreeT::Ptr newTree(
new DistTreeT(exBandWidth));
1901 newTree->merge(mDistGrid->tree());
1902 mDistGrid->setTree(newTree);
1908 typedef internal::OffsetOp<DistValueT> OffsetOp;
1909 typedef internal::RenormOp<DistGridT, DistValueT> RenormOp;
1910 typedef internal::MinOp<DistTreeT, DistValueT> MinOp;
1911 typedef internal::MergeBufferOp<DistTreeT, DistValueT> MergeBufferOp;
1913 tree::LeafManager<DistTreeT> leafs(mDistGrid->tree(), 1);
1915 const DistValueT offset = 0.8 * voxelSize;
1917 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1919 OffsetOp offsetOp(-offset);
1920 LeafTransformer<DistTreeT, OffsetOp> offsetXform(leafs, offsetOp);
1921 offsetXform.runParallel();
1923 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1925 RenormOp renormOp(*mDistGrid, leafs, voxelSize);
1926 LeafTransformer<DistTreeT, RenormOp> renormXform(leafs, renormOp);
1927 renormXform.runParallel();
1930 LeafTransformer<DistTreeT, MinOp> minXform(leafs, minOp);
1931 minXform.runParallel();
1933 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1935 offsetOp.resetOffset(offset - internal::Tolerance<DistValueT>::epsilon());
1936 offsetXform.runParallel();
1939 mIntersectingVoxelsGrid->clear();
1942 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1945 const DistValueT minWidth = voxelSize * 2.0;
1946 if (inBandWidth > minWidth || exBandWidth > minWidth) {
1949 StencilTreeT maskTree(
false);
1950 tree::ValueAccessor<StencilTreeT> acc(maskTree);
1951 maskTree.topologyUnion(mDistGrid->tree());
1955 typedef typename DistTreeT::LeafNodeType DistLeafType;
1957 std::vector<DistLeafType*> distLeafs;
1958 distLeafs.reserve(mDistGrid->tree().leafCount());
1960 typename DistTreeT::LeafIter iter = mDistGrid->tree().beginLeaf();
1961 for ( ; iter; ++iter) distLeafs.push_back(iter.getLeaf());
1963 tree::ValueAccessor<DistTreeT> distAcc(mDistGrid->tree());
1965 DistValueT leafSize = DistValueT(DistLeafType::DIM - 1) * voxelSize;
1967 const double inLeafsRatio = double(inBandWidth) / double(leafSize);
1969 if (
double(inLeafs) > (inLeafsRatio + 1.0)) {
1970 inLeafs = size_t(std::ceil(inLeafsRatio)) + 1;
1972 size_t exLeafs = size_t(std::ceil(exBandWidth / leafSize)) + 1;
1973 size_t numLeafs =
std::max(inLeafs, exLeafs);
1975 for (
size_t i = 0; i < numLeafs; ++i) {
1977 if (mInterrupter && mInterrupter->wasInterrupted())
return;
1979 std::vector<DistLeafType*> newDistLeafs;
1980 newDistLeafs.reserve(2 * distLeafs.size());
1982 for (
size_t n = 0, N = distLeafs.size(); n < N; ++n) {
1984 Coord ijk = distLeafs[n]->getOrigin();
1986 const bool inside = distLeafs[n]->getValue(ijk) < DistValueT(0.0);
1988 if (inside && !(i < inLeafs))
continue;
1989 else if (!inside && !(i < exLeafs))
continue;
1992 if (distAcc.probeLeaf(ijk) == NULL) {
1993 newDistLeafs.push_back(distAcc.touchLeaf(ijk));
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));
2009 ijk[0] += DistLeafType::DIM;
2010 if (distAcc.probeLeaf(ijk) == NULL) {
2011 newDistLeafs.push_back(distAcc.touchLeaf(ijk));
2014 ijk[0] -= DistLeafType::DIM;
2015 ijk[1] += DistLeafType::DIM;
2016 if (distAcc.probeLeaf(ijk) == NULL) {
2017 newDistLeafs.push_back(distAcc.touchLeaf(ijk));
2020 ijk[1] -= DistLeafType::DIM;
2021 ijk[2] += DistLeafType::DIM;
2022 if (distAcc.probeLeaf(ijk) == NULL) {
2023 newDistLeafs.push_back(distAcc.touchLeaf(ijk));
2027 if (newDistLeafs.empty())
break;
2028 distLeafs.swap(newDistLeafs);
2032 if (mInterrupter && mInterrupter->wasInterrupted())
return;
2034 mIndexGrid->tree().topologyUnion(mDistGrid->tree());
2036 while (maskTree.activeVoxelCount() > 0) {
2038 if (mInterrupter && mInterrupter->wasInterrupted())
break;
2041 tree::LeafManager<StencilTreeT> leafs(maskTree);
2043 internal::ExpandNB<DistTreeT> expand(pointList, polygonList, mDistGrid->tree(),
2044 mIndexGrid->tree(), maskTree, leafs, exBandWidth, inBandWidth, voxelSize);
2046 expand.runParallel();
2052 const DistValueT minTrimWidth = voxelSize * 4.0;
2053 if (inBandWidth < minTrimWidth || exBandWidth < minTrimWidth) {
2059 tree::LeafManager<DistTreeT> leafs(mDistGrid->tree());
2061 typedef internal::TrimOp<DistValueT> TrimOp;
2063 TrimOp op(exBandWidth, unsignedDistField ? exBandWidth : inBandWidth);
2064 LeafTransformer<DistTreeT, TrimOp> transform(leafs, op);
2065 transform.runParallel();
2068 if (mInterrupter && mInterrupter->wasInterrupted())
return;
2070 mDistGrid->tree().pruneLevelSet();
2077 #endif // OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED