44 #ifndef OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
45 #define OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
47 #include <openvdb/Types.h>
48 #include <openvdb/math/FiniteDifference.h>
49 #include <openvdb/math/Proximity.h>
50 #include <openvdb/util/NullInterrupter.h>
51 #include <openvdb/util/Util.h>
57 #include <tbb/blocked_range.h>
58 #include <tbb/enumerable_thread_specific.h>
59 #include <tbb/parallel_for.h>
60 #include <tbb/parallel_reduce.h>
61 #include <tbb/partitioner.h>
62 #include <tbb/task_group.h>
63 #include <tbb/task_scheduler_init.h>
64 #include <tbb/tick_count.h>
66 #include <boost/integer_traits.hpp>
67 #include <boost/math/special_functions/fpclassify.hpp>
68 #include <boost/scoped_array.hpp>
137 template <
typename Gr
idType,
typename MeshDataAdapter>
138 inline typename GridType::Ptr
142 float exteriorBandWidth = 3.0f,
143 float interiorBandWidth = 3.0f,
163 template <
typename Gr
idType,
typename MeshDataAdapter,
typename Interrupter>
164 inline typename GridType::Ptr
166 Interrupter& interrupter,
169 float exteriorBandWidth = 3.0f,
170 float interiorBandWidth = 3.0f,
188 template<
typename Po
intType,
typename PolygonType>
192 const std::vector<PolygonType>& polygons)
193 : mPointArray(points.empty() ? NULL : &points[0])
194 , mPointArraySize(points.size())
195 , mPolygonArray(polygons.empty() ? NULL : &polygons[0])
196 , mPolygonArraySize(polygons.size())
201 const PolygonType* polygonArray,
size_t polygonArraySize)
202 : mPointArray(pointArray)
203 , mPointArraySize(pointArraySize)
204 , mPolygonArray(polygonArray)
205 , mPolygonArraySize(polygonArraySize)
214 return (PolygonType::size == 3 || mPolygonArray[n][3] ==
util::INVALID_IDX) ? 3 : 4;
220 const PointType& p = mPointArray[mPolygonArray[n][int(v)]];
221 pos[0] = double(p[0]);
222 pos[1] = double(p[1]);
223 pos[2] = double(p[2]);
227 PointType
const *
const mPointArray;
228 size_t const mPointArraySize;
229 PolygonType
const *
const mPolygonArray;
230 size_t const mPolygonArraySize;
255 template<
typename Gr
idType>
256 inline typename GridType::Ptr
258 const openvdb::math::Transform& xform,
259 const std::vector<Vec3s>& points,
260 const std::vector<Vec3I>& triangles,
279 template<
typename Gr
idType>
280 inline typename GridType::Ptr
282 const openvdb::math::Transform& xform,
283 const std::vector<Vec3s>& points,
284 const std::vector<Vec4I>& quads,
304 template<
typename Gr
idType>
305 inline typename GridType::Ptr
307 const openvdb::math::Transform& xform,
308 const std::vector<Vec3s>& points,
309 const std::vector<Vec3I>& triangles,
310 const std::vector<Vec4I>& quads,
332 template<
typename Gr
idType>
333 inline typename GridType::Ptr
335 const openvdb::math::Transform& xform,
336 const std::vector<Vec3s>& points,
337 const std::vector<Vec3I>& triangles,
338 const std::vector<Vec4I>& quads,
357 template<
typename Gr
idType>
358 inline typename GridType::Ptr
360 const openvdb::math::Transform& xform,
361 const std::vector<Vec3s>& points,
362 const std::vector<Vec3I>& triangles,
363 const std::vector<Vec4I>& quads,
376 template<
typename Gr
idType,
typename VecType>
377 inline typename GridType::Ptr
379 const openvdb::math::Transform& xform,
392 template <
typename FloatTreeT>
410 : mXDist(dist), mYDist(dist), mZDist(dist)
453 void convert(
const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList);
458 void getEdgeData(Accessor& acc,
const Coord& ijk,
459 std::vector<Vec3d>& points, std::vector<Index32>& primitives);
478 namespace mesh_to_volume_internal {
480 template<
typename Po
intType>
485 : mPointsIn(pointsIn), mPointsOut(pointsOut), mXform(&xform)
489 void operator()(
const tbb::blocked_range<size_t>& range)
const {
493 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
495 const PointType& wsP = mPointsIn[n];
496 pos[0] = double(wsP[0]);
497 pos[1] = double(wsP[1]);
498 pos[2] = double(wsP[2]);
500 pos = mXform->worldToIndex(pos);
502 PointType& isP = mPointsOut[n];
503 isP[0] =
typename PointType::value_type(pos[0]);
504 isP[1] =
typename PointType::value_type(pos[1]);
505 isP[2] =
typename PointType::value_type(pos[2]);
515 template<
typename ValueType>
518 static ValueType
epsilon() {
return ValueType(1e-7); }
526 template<
typename TreeType>
537 LeafNodeType ** rhsDistNodes, Int32LeafNodeType ** rhsIdxNodes)
538 : mDistTree(&lhsDistTree)
539 , mIdxTree(&lhsIdxTree)
540 , mRhsDistNodes(rhsDistNodes)
541 , mRhsIdxNodes(rhsIdxNodes)
545 void operator()(
const tbb::blocked_range<size_t>& range)
const {
550 typedef typename LeafNodeType::ValueType DistValueType;
551 typedef typename Int32LeafNodeType::ValueType IndexValueType;
553 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
555 const Coord& origin = mRhsDistNodes[n]->origin();
557 LeafNodeType* lhsDistNode = distAcc.
probeLeaf(origin);
558 Int32LeafNodeType* lhsIdxNode = idxAcc.
probeLeaf(origin);
560 DistValueType* lhsDistData = lhsDistNode->buffer().data();
561 IndexValueType* lhsIdxData = lhsIdxNode->buffer().data();
563 const DistValueType* rhsDistData = mRhsDistNodes[n]->buffer().data();
564 const IndexValueType* rhsIdxData = mRhsIdxNodes[n]->buffer().data();
567 for (
Index32 offset = 0; offset < LeafNodeType::SIZE; ++offset) {
571 const DistValueType& lhsValue = lhsDistData[offset];
572 const DistValueType& rhsValue = rhsDistData[offset];
574 if (rhsValue < lhsValue) {
575 lhsDistNode->setValueOn(offset, rhsValue);
576 lhsIdxNode->setValueOn(offset, rhsIdxData[offset]);
578 lhsIdxNode->setValueOn(offset,
579 std::min(lhsIdxData[offset], rhsIdxData[offset]));
584 delete mRhsDistNodes[n];
585 delete mRhsIdxNodes[n];
591 TreeType *
const mDistTree;
592 Int32TreeType *
const mIdxTree;
594 LeafNodeType **
const mRhsDistNodes;
595 Int32LeafNodeType **
const mRhsIdxNodes;
602 template<
typename TreeType>
608 : mNodes(nodes.empty() ? NULL : &nodes[0]), mCoordinates(coordinates)
612 void operator()(
const tbb::blocked_range<size_t>& range)
const {
613 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
614 Coord& origin =
const_cast<Coord&
>(mNodes[n]->origin());
615 mCoordinates[n] = origin;
616 origin[0] =
static_cast<int>(n);
625 template<
typename TreeType>
631 : mNodes(nodes.empty() ? NULL : &nodes[0]), mCoordinates(coordinates)
635 void operator()(
const tbb::blocked_range<size_t>& range)
const {
636 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
637 Coord& origin =
const_cast<Coord&
>(mNodes[n]->origin());
638 origin[0] = mCoordinates[n][0];
647 template<
typename TreeType>
654 size_t* offsets,
size_t numNodes,
const CoordBBox& bbox)
656 , mCoordinates(coordinates)
658 , mNumNodes(numNodes)
663 void operator()(
const tbb::blocked_range<size_t>& range)
const {
665 size_t* offsetsNextX = mOffsets;
666 size_t* offsetsPrevX = mOffsets + mNumNodes;
667 size_t* offsetsNextY = mOffsets + mNumNodes * 2;
668 size_t* offsetsPrevY = mOffsets + mNumNodes * 3;
669 size_t* offsetsNextZ = mOffsets + mNumNodes * 4;
670 size_t* offsetsPrevZ = mOffsets + mNumNodes * 5;
675 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
676 const Coord& origin = mCoordinates[n];
677 offsetsNextX[n] = findNeighbourNode(acc, origin, Coord(LeafNodeType::DIM, 0, 0));
678 offsetsPrevX[n] = findNeighbourNode(acc, origin, Coord(-LeafNodeType::DIM, 0, 0));
679 offsetsNextY[n] = findNeighbourNode(acc, origin, Coord(0, LeafNodeType::DIM, 0));
680 offsetsPrevY[n] = findNeighbourNode(acc, origin, Coord(0, -LeafNodeType::DIM, 0));
681 offsetsNextZ[n] = findNeighbourNode(acc, origin, Coord(0, 0, LeafNodeType::DIM));
682 offsetsPrevZ[n] = findNeighbourNode(acc, origin, Coord(0, 0, -LeafNodeType::DIM));
688 Coord ijk = start + step;
689 CoordBBox bbox(mBBox);
691 while (bbox.isInside(ijk)) {
693 if (node)
return static_cast<size_t>(node->origin()[0]);
697 return boost::integer_traits<size_t>::const_max;
705 TreeType
const *
const mTree;
706 Coord
const *
const mCoordinates;
707 size_t *
const mOffsets;
709 const size_t mNumNodes;
710 const CoordBBox mBBox;
714 template<
typename TreeType>
717 enum { INVALID_OFFSET = boost::integer_traits<size_t>::const_max };
725 mLeafNodes.reserve(tree.leafCount());
726 tree.getNodes(mLeafNodes);
728 if (mLeafNodes.empty())
return;
731 tree.evalLeafBoundingBox(bbox);
733 const tbb::blocked_range<size_t> range(0, mLeafNodes.size());
737 boost::scoped_array<Coord> coordinates(
new Coord[mLeafNodes.size()]);
741 mOffsets.reset(
new size_t[mLeafNodes.size() * 6]);
744 tbb::parallel_for(range,
751 size_t size()
const {
return mLeafNodes.size(); }
753 std::vector<LeafNodeType*>&
nodes() {
return mLeafNodes; }
754 const std::vector<LeafNodeType*>&
nodes()
const {
return mLeafNodes; }
758 const size_t*
offsetsPrevX()
const {
return mOffsets.get() + mLeafNodes.size(); }
760 const size_t*
offsetsNextY()
const {
return mOffsets.get() + mLeafNodes.size() * 2; }
761 const size_t*
offsetsPrevY()
const {
return mOffsets.get() + mLeafNodes.size() * 3; }
763 const size_t*
offsetsNextZ()
const {
return mOffsets.get() + mLeafNodes.size() * 4; }
764 const size_t*
offsetsPrevZ()
const {
return mOffsets.get() + mLeafNodes.size() * 5; }
767 std::vector<LeafNodeType*> mLeafNodes;
768 boost::scoped_array<size_t> mOffsets;
772 template<
typename TreeType>
784 : mStartNodeIndices(startNodeIndices.empty() ? NULL : &startNodeIndices[0])
785 , mConnectivity(&connectivity)
790 void operator()(
const tbb::blocked_range<size_t>& range)
const {
792 std::vector<LeafNodeType*>& nodes = mConnectivity->nodes();
795 size_t idxA = 0, idxB = 1;
798 const size_t* nextOffsets = mConnectivity->offsetsNextZ();
799 const size_t* prevOffsets = mConnectivity->offsetsPrevZ();
805 step = LeafNodeType::DIM;
807 nextOffsets = mConnectivity->offsetsNextY();
808 prevOffsets = mConnectivity->offsetsPrevY();
810 }
else if (mAxis ==
X_AXIS) {
814 step = LeafNodeType::DIM * LeafNodeType::DIM;
816 nextOffsets = mConnectivity->offsetsNextX();
817 prevOffsets = mConnectivity->offsetsPrevX();
825 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
827 size_t startOffset = mStartNodeIndices[n];
828 size_t lastOffset = startOffset;
832 for (a = 0; a < int(LeafNodeType::DIM); ++a) {
833 for (b = 0; b < int(LeafNodeType::DIM); ++b) {
835 pos = LeafNodeType::coordToOffset(ijk);
836 size_t offset = startOffset;
839 while ( offset != ConnectivityTable::INVALID_OFFSET &&
840 traceVoxelLine(*nodes[offset], pos, step) ) {
843 offset = nextOffsets[offset];
848 while (offset != ConnectivityTable::INVALID_OFFSET) {
850 offset = nextOffsets[offset];
855 pos += step * (LeafNodeType::DIM - 1);
856 while ( offset != ConnectivityTable::INVALID_OFFSET &&
857 traceVoxelLine(*nodes[offset], pos, -step)) {
858 offset = prevOffsets[offset];
868 ValueType* data = node.buffer().data();
870 bool isOutside =
true;
872 for (
Index i = 0; i < LeafNodeType::DIM; ++i) {
874 ValueType& dist = data[pos];
876 if (dist < ValueType(0.0)) {
880 if (!(dist > ValueType(0.75))) isOutside =
false;
882 if (isOutside) dist = ValueType(-dist);
893 size_t const *
const mStartNodeIndices;
894 ConnectivityTable *
const mConnectivity;
900 template<
typename LeafNodeType>
904 typedef typename LeafNodeType::ValueType ValueType;
905 typedef std::deque<Index> Queue;
908 ValueType* data = node.buffer().data();
912 for (
Index pos = 0; pos < LeafNodeType::SIZE; ++pos) {
913 if (data[pos] < 0.0) seedPoints.push_back(pos);
916 if (seedPoints.empty())
return;
919 for (Queue::iterator it = seedPoints.begin(); it != seedPoints.end(); ++it) {
920 ValueType& dist = data[*it];
927 Index pos(0), nextPos(0);
929 while (!seedPoints.empty()) {
931 pos = seedPoints.back();
932 seedPoints.pop_back();
934 ValueType& dist = data[pos];
936 if (!(dist < ValueType(0.0))) {
940 ijk = LeafNodeType::offsetToLocalCoord(pos);
943 nextPos = pos - LeafNodeType::DIM * LeafNodeType::DIM;
944 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
947 if (ijk[0] != (LeafNodeType::DIM - 1)) {
948 nextPos = pos + LeafNodeType::DIM * LeafNodeType::DIM;
949 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
953 nextPos = pos - LeafNodeType::DIM;
954 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
957 if (ijk[1] != (LeafNodeType::DIM - 1)) {
958 nextPos = pos + LeafNodeType::DIM;
959 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
964 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
967 if (ijk[2] != (LeafNodeType::DIM - 1)) {
969 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
976 template<
typename LeafNodeType>
980 bool updatedNode =
false;
982 typedef typename LeafNodeType::ValueType ValueType;
983 ValueType* data = node.buffer().data();
987 bool updatedSign =
true;
988 while (updatedSign) {
992 for (
Index pos = 0; pos < LeafNodeType::SIZE; ++pos) {
994 ValueType& dist = data[pos];
996 if (!(dist < ValueType(0.0)) && dist > ValueType(0.75)) {
998 ijk = LeafNodeType::offsetToLocalCoord(pos);
1001 if (ijk[2] != 0 && data[pos - 1] < ValueType(0.0)) {
1003 dist = ValueType(-dist);
1006 }
else if (ijk[2] != (LeafNodeType::DIM - 1) && data[pos + 1] < ValueType(0.0)) {
1008 dist = ValueType(-dist);
1011 }
else if (ijk[1] != 0 && data[pos - LeafNodeType::DIM] < ValueType(0.0)) {
1013 dist = ValueType(-dist);
1016 }
else if (ijk[1] != (LeafNodeType::DIM - 1) && data[pos + LeafNodeType::DIM] < ValueType(0.0)) {
1018 dist = ValueType(-dist);
1021 }
else if (ijk[0] != 0 && data[pos - LeafNodeType::DIM * LeafNodeType::DIM] < ValueType(0.0)) {
1023 dist = ValueType(-dist);
1026 }
else if (ijk[0] != (LeafNodeType::DIM - 1) && data[pos + LeafNodeType::DIM * LeafNodeType::DIM] < ValueType(0.0)) {
1028 dist = ValueType(-dist);
1033 updatedNode |= updatedSign;
1040 template<
typename TreeType>
1048 : mNodes(nodes.empty() ? NULL : &nodes[0])
1049 , mChangedNodeMask(changedNodeMask)
1054 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1055 if (mChangedNodeMask[n]) {
1057 mChangedNodeMask[n] =
scanFill(*mNodes[n]);
1067 template<
typename ValueType>
1070 FillArray(ValueType* array,
const ValueType v) : mArray(array), mValue(v) { }
1073 const ValueType v = mValue;
1074 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1084 template<
typename ValueType>
1086 fillArray(ValueType* array,
const ValueType val,
const size_t length)
1088 const size_t grainSize = length / tbb::task_scheduler_init::default_num_threads();
1089 const tbb::blocked_range<size_t> range(0, length, grainSize);
1094 template<
typename TreeType>
1101 SyncVoxelMask(std::vector<LeafNodeType*>& nodes,
const bool* changedNodeMask,
bool* changedVoxelMask)
1102 : mNodes(nodes.empty() ? NULL : &nodes[0])
1103 , mChangedNodeMask(changedNodeMask)
1104 , mChangedVoxelMask(changedVoxelMask)
1109 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1111 if (mChangedNodeMask[n]) {
1112 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1114 ValueType* data = mNodes[n]->buffer().data();
1116 for (
Index pos = 0; pos < LeafNodeType::SIZE; ++pos) {
1118 data[pos] = ValueType(-data[pos]);
1132 template<
typename TreeType>
1140 SeedPoints(ConnectivityTable& connectivity,
bool* changedNodeMask,
bool* nodeMask,
bool* changedVoxelMask)
1141 : mConnectivity(&connectivity)
1142 , mChangedNodeMask(changedNodeMask)
1143 , mNodeMask(nodeMask)
1144 , mChangedVoxelMask(changedVoxelMask)
1150 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1152 if (!mChangedNodeMask[n]) {
1154 bool changedValue =
false;
1156 changedValue |= processZ(n,
true);
1157 changedValue |= processZ(n,
false);
1159 changedValue |= processY(n,
true);
1160 changedValue |= processY(n,
false);
1162 changedValue |= processX(n,
true);
1163 changedValue |= processX(n,
false);
1165 mNodeMask[n] = changedValue;
1173 const size_t offset = firstFace ? mConnectivity->offsetsPrevZ()[n] : mConnectivity->offsetsNextZ()[n];
1174 if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1176 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1178 const ValueType* lhsData = mConnectivity->nodes()[n]->buffer().data();
1179 const ValueType* rhsData = mConnectivity->nodes()[offset]->buffer().data();
1181 const Index lastOffset = LeafNodeType::DIM - 1;
1182 const Index lhsOffset = firstFace ? 0 : lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1184 Index tmpPos(0), pos(0);
1185 bool changedValue =
false;
1187 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
1188 tmpPos = x << (2 * LeafNodeType::LOG2DIM);
1189 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
1190 pos = tmpPos + (y << LeafNodeType::LOG2DIM);
1192 if (lhsData[pos + lhsOffset] > ValueType(0.75)) {
1193 if (rhsData[pos + rhsOffset] < ValueType(0.0)) {
1194 changedValue =
true;
1195 mask[pos + lhsOffset] =
true;
1201 return changedValue;
1209 const size_t offset = firstFace ? mConnectivity->offsetsPrevY()[n] : mConnectivity->offsetsNextY()[n];
1210 if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1212 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1214 const ValueType* lhsData = mConnectivity->nodes()[n]->buffer().data();
1215 const ValueType* rhsData = mConnectivity->nodes()[offset]->buffer().data();
1217 const Index lastOffset = LeafNodeType::DIM * (LeafNodeType::DIM - 1);
1218 const Index lhsOffset = firstFace ? 0 : lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1220 Index tmpPos(0), pos(0);
1221 bool changedValue =
false;
1223 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
1224 tmpPos = x << (2 * LeafNodeType::LOG2DIM);
1225 for (
Index z = 0; z < LeafNodeType::DIM; ++z) {
1228 if (lhsData[pos + lhsOffset] > ValueType(0.75)) {
1229 if (rhsData[pos + rhsOffset] < ValueType(0.0)) {
1230 changedValue =
true;
1231 mask[pos + lhsOffset] =
true;
1237 return changedValue;
1245 const size_t offset = firstFace ? mConnectivity->offsetsPrevX()[n] : mConnectivity->offsetsNextX()[n];
1246 if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1248 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1250 const ValueType* lhsData = mConnectivity->nodes()[n]->buffer().data();
1251 const ValueType* rhsData = mConnectivity->nodes()[offset]->buffer().data();
1253 const Index lastOffset = LeafNodeType::DIM * LeafNodeType::DIM * (LeafNodeType::DIM - 1);
1254 const Index lhsOffset = firstFace ? 0 : lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1256 Index tmpPos(0), pos(0);
1257 bool changedValue =
false;
1259 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
1260 tmpPos = y << LeafNodeType::LOG2DIM;
1261 for (
Index z = 0; z < LeafNodeType::DIM; ++z) {
1264 if (lhsData[pos + lhsOffset] > ValueType(0.75)) {
1265 if (rhsData[pos + rhsOffset] < ValueType(0.0)) {
1266 changedValue =
true;
1267 mask[pos + lhsOffset] =
true;
1273 return changedValue;
1288 template<
typename TreeType,
typename MeshDataAdapter>
1296 typedef std::pair<boost::shared_array<Vec3d>, boost::shared_array<bool> >
LocalData;
1300 std::vector<LeafNodeType*>& distNodes,
1301 const TreeType& distTree,
1302 const Int32TreeType& indexTree,
1304 : mDistNodes(distNodes.empty() ? NULL : &distNodes[0])
1305 , mDistTree(&distTree)
1306 , mIndexTree(&indexTree)
1308 , mLocalDataTable(new LocalDataTable())
1320 Index xPos(0), yPos(0);
1321 Coord ijk, nijk, nodeMin, nodeMax;
1322 Vec3d cp, xyz, nxyz, dir1, dir2;
1324 LocalData& localData = mLocalDataTable->local();
1326 boost::shared_array<Vec3d>& points = localData.first;
1327 if (!points) points.reset(
new Vec3d[LeafNodeType::SIZE * 2]);
1329 boost::shared_array<bool>& mask = localData.second;
1330 if (!mask) mask.reset(
new bool[LeafNodeType::SIZE]);
1333 typename LeafNodeType::ValueOnCIter it;
1335 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1337 LeafNodeType& node = *mDistNodes[n];
1338 ValueType* data = node.buffer().data();
1340 const Int32LeafNodeType* idxNode = idxAcc.
probeConstLeaf(node.origin());
1341 const Int32* idxData = idxNode->buffer().data();
1343 nodeMin = node.origin();
1344 nodeMax = nodeMin.offsetBy(LeafNodeType::DIM - 1);
1347 memset(mask.get(), 0,
sizeof(bool) * LeafNodeType::SIZE);
1349 for (it = node.cbeginValueOn(); it; ++it) {
1350 Index pos = it.pos();
1352 ValueType& dist = data[pos];
1353 if (dist < 0.0 || dist > 0.75)
continue;
1355 ijk = node.offsetToGlobalCoord(pos);
1357 xyz[0] = double(ijk[0]);
1358 xyz[1] = double(ijk[1]);
1359 xyz[2] = double(ijk[2]);
1365 bool flipSign =
false;
1367 for (nijk[0] = bbox.min()[0]; nijk[0] <= bbox.max()[0] && !flipSign; ++nijk[0]) {
1368 xPos = (nijk[0] & (LeafNodeType::DIM - 1u)) << (2 * LeafNodeType::LOG2DIM);
1369 for (nijk[1] = bbox.min()[1]; nijk[1] <= bbox.max()[1] && !flipSign; ++nijk[1]) {
1370 yPos = xPos + ((nijk[1] & (LeafNodeType::DIM - 1u)) << LeafNodeType::LOG2DIM);
1371 for (nijk[2] = bbox.min()[2]; nijk[2] <= bbox.max()[2]; ++nijk[2]) {
1372 pos = yPos + (nijk[2] & (LeafNodeType::DIM - 1u));
1374 const Int32& polyIdx = idxData[pos];
1378 const Index pointIndex = pos * 2;
1384 nxyz[0] = double(nijk[0]);
1385 nxyz[1] = double(nijk[1]);
1386 nxyz[2] = double(nijk[2]);
1388 Vec3d& point = points[pointIndex];
1390 point = closestPoint(nxyz, polyIdx);
1392 Vec3d& direction = points[pointIndex + 1];
1393 direction = nxyz - point;
1397 dir1 = xyz - points[pointIndex];
1400 if (points[pointIndex + 1].dot(dir1) > 0.0) {
1411 for (
Int32 m = 0; m < 26; ++m) {
1414 if (!bbox.isInside(nijk) && distAcc.
probeValue(nijk, nval) && nval < -0.75) {
1415 nxyz[0] = double(nijk[0]);
1416 nxyz[1] = double(nijk[1]);
1417 nxyz[2] = double(nijk[2]);
1419 cp = closestPoint(nxyz, idxAcc.
getValue(nijk));
1427 if (dir2.dot(dir1) > 0.0) {
1443 Vec3d a, b, c, cp, uvw;
1445 const size_t polygon = size_t(polyIdx);
1446 mMesh->getIndexSpacePoint(polygon, 0, a);
1447 mMesh->getIndexSpacePoint(polygon, 1, b);
1448 mMesh->getIndexSpacePoint(polygon, 2, c);
1452 if (4 == mMesh->vertexCount(polygon)) {
1454 mMesh->getIndexSpacePoint(polygon, 3, b);
1458 if ((center - c).lengthSqr() < (center - cp).lengthSqr()) {
1467 LeafNodeType **
const mDistNodes;
1468 TreeType
const *
const mDistTree;
1469 Int32TreeType
const *
const mIndexTree;
1472 boost::shared_ptr<LocalDataTable> mLocalDataTable;
1479 template<
typename LeafNodeType>
1483 typedef LeafNodeType NodeT;
1485 const Coord ijk = NodeT::offsetToLocalCoord(pos);
1489 mask[0] = ijk[0] != (NodeT::DIM - 1);
1491 mask[1] = ijk[0] != 0;
1493 mask[2] = ijk[1] != (NodeT::DIM - 1);
1495 mask[3] = ijk[1] != 0;
1497 mask[4] = ijk[2] != (NodeT::DIM - 1);
1499 mask[5] = ijk[2] != 0;
1503 mask[6] = mask[0] && mask[5];
1505 mask[7] = mask[1] && mask[5];
1507 mask[8] = mask[0] && mask[4];
1509 mask[9] = mask[1] && mask[4];
1511 mask[10] = mask[0] && mask[2];
1513 mask[11] = mask[1] && mask[2];
1515 mask[12] = mask[0] && mask[3];
1517 mask[13] = mask[1] && mask[3];
1519 mask[14] = mask[3] && mask[4];
1521 mask[15] = mask[3] && mask[5];
1523 mask[16] = mask[2] && mask[4];
1525 mask[17] = mask[2] && mask[5];
1529 mask[18] = mask[1] && mask[3] && mask[5];
1531 mask[19] = mask[1] && mask[3] && mask[4];
1533 mask[20] = mask[0] && mask[3] && mask[4];
1535 mask[21] = mask[0] && mask[3] && mask[5];
1537 mask[22] = mask[1] && mask[2] && mask[5];
1539 mask[23] = mask[1] && mask[2] && mask[4];
1541 mask[24] = mask[0] && mask[2] && mask[4];
1543 mask[25] = mask[0] && mask[2] && mask[5];
1547 template<
typename Compare,
typename LeafNodeType>
1551 typedef LeafNodeType NodeT;
1554 if (mask[5] && Compare::check(data[pos - 1]))
return true;
1556 if (mask[4] && Compare::check(data[pos + 1]))
return true;
1558 if (mask[3] && Compare::check(data[pos - NodeT::DIM]))
return true;
1560 if (mask[2] && Compare::check(data[pos + NodeT::DIM]))
return true;
1562 if (mask[1] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM]))
return true;
1564 if (mask[0] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM]))
return true;
1566 if (mask[6] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM]))
return true;
1568 if (mask[7] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - 1]))
return true;
1570 if (mask[8] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + 1]))
return true;
1572 if (mask[9] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + 1]))
return true;
1574 if (mask[10] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM]))
return true;
1576 if (mask[11] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM]))
return true;
1578 if (mask[12] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM]))
return true;
1580 if (mask[13] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM]))
return true;
1582 if (mask[14] && Compare::check(data[pos - NodeT::DIM + 1]))
return true;
1584 if (mask[15] && Compare::check(data[pos - NodeT::DIM - 1]))
return true;
1586 if (mask[16] && Compare::check(data[pos + NodeT::DIM + 1]))
return true;
1588 if (mask[17] && Compare::check(data[pos + NodeT::DIM - 1]))
return true;
1590 if (mask[18] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM - 1]))
return true;
1592 if (mask[19] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM + 1]))
return true;
1594 if (mask[20] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM + 1]))
return true;
1596 if (mask[21] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM - 1]))
return true;
1598 if (mask[22] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM - 1]))
return true;
1600 if (mask[23] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM + 1]))
return true;
1602 if (mask[24] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM + 1]))
return true;
1604 if (mask[25] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM - 1]))
return true;
1610 template<
typename Compare,
typename AccessorType>
1614 for (
Int32 m = 0; m < 26; ++m) {
1624 template<
typename TreeType>
1630 struct IsNegative {
static bool check(
const ValueType v) {
return v < ValueType(0.0); } };
1634 , mNodes(nodes.empty() ? NULL : &nodes[0])
1641 bool neighbourMask[26];
1643 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1645 LeafNodeType& node = *mNodes[n];
1646 ValueType* data = node.buffer().data();
1648 typename LeafNodeType::ValueOnCIter it;
1649 for (it = node.cbeginValueOn(); it; ++it) {
1651 const Index pos = it.pos();
1653 ValueType& dist = data[pos];
1654 if (dist < 0.0 || dist > 0.75)
continue;
1657 maskNodeInternalNeighbours<LeafNodeType>(pos, neighbourMask);
1659 const bool hasNegativeNeighbour =
1660 checkNeighbours<IsNegative, LeafNodeType>(pos, data, neighbourMask) ||
1661 checkNeighbours<IsNegative>(node.offsetToGlobalCoord(pos), acc, neighbourMask);
1663 if (!hasNegativeNeighbour) {
1676 template<
typename TreeType>
1683 struct Comp {
static bool check(
const ValueType v) {
return !(v > ValueType(0.75)); } };
1686 TreeType& distTree, Int32TreeType& indexTree)
1687 : mNodes(nodes.empty() ? NULL : &nodes[0])
1688 , mDistTree(&distTree)
1689 , mIndexTree(&indexTree)
1697 bool neighbourMask[26];
1699 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1701 LeafNodeType& distNode = *mNodes[n];
1702 ValueType* data = distNode.buffer().data();
1704 typename Int32TreeType::LeafNodeType* idxNode =
1707 typename LeafNodeType::ValueOnCIter it;
1708 for (it = distNode.cbeginValueOn(); it; ++it) {
1710 const Index pos = it.pos();
1712 if (!(data[pos] > 0.75))
continue;
1715 maskNodeInternalNeighbours<LeafNodeType>(pos, neighbourMask);
1717 const bool hasBoundaryNeighbour =
1718 checkNeighbours<Comp, LeafNodeType>(pos, data, neighbourMask) ||
1719 checkNeighbours<Comp>(distNode.offsetToGlobalCoord(pos), distAcc, neighbourMask);
1721 if (!hasBoundaryNeighbour) {
1722 distNode.setValueOff(pos);
1723 idxNode->setValueOff(pos);
1738 template<
typename NodeType>
1745 typedef typename NodeType::NodeMaskType NodeMaskType;
1747 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1748 const_cast<NodeMaskType&
>(mNodes[n]->getChildMask()).setOff();
1756 template<
typename TreeType>
1760 typedef typename TreeType::RootNodeType RootNodeType;
1761 typedef typename RootNodeType::NodeChainType NodeChainType;
1762 typedef typename boost::mpl::at<NodeChainType, boost::mpl::int_<1> >::type InternalNodeType;
1764 std::vector<InternalNodeType*> nodes;
1765 tree.getNodes(nodes);
1767 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
1772 template<
typename TreeType>
1778 std::vector<LeafNodeType*>& overlappingNodes)
1779 : mLhsTree(&lhsTree)
1780 , mRhsTree(&rhsTree)
1781 , mNodes(&overlappingNodes)
1787 std::vector<LeafNodeType*> rhsLeafNodes;
1789 rhsLeafNodes.reserve(mRhsTree->leafCount());
1792 mRhsTree->stealNodes(rhsLeafNodes);
1796 for (
size_t n = 0, N = rhsLeafNodes.size(); n < N; ++n) {
1797 if (!acc.
probeLeaf(rhsLeafNodes[n]->origin())) {
1800 mNodes->push_back(rhsLeafNodes[n]);
1806 TreeType *
const mLhsTree;
1807 TreeType *
const mRhsTree;
1808 std::vector<LeafNodeType*> *
const mNodes;
1812 template<
typename DistTreeType,
typename IndexTreeType>
1815 DistTreeType& rhsDist, IndexTreeType& rhsIdx)
1817 typedef typename DistTreeType::LeafNodeType DistLeafNodeType;
1818 typedef typename IndexTreeType::LeafNodeType IndexLeafNodeType;
1820 std::vector<DistLeafNodeType*> overlappingDistNodes;
1821 std::vector<IndexLeafNodeType*> overlappingIdxNodes;
1824 tbb::task_group tasks;
1830 if (!overlappingDistNodes.empty() && !overlappingIdxNodes.empty()) {
1831 tbb::parallel_for(tbb::blocked_range<size_t>(0, overlappingDistNodes.size()),
1843 template<
typename TreeType>
1846 typedef boost::scoped_ptr<VoxelizationData>
Ptr;
1850 typedef typename TreeType::template ValueConverter<unsigned char>::Type
UCharTreeType;
1858 : distTree(std::numeric_limits<ValueType>::
max())
1861 , indexAcc(indexTree)
1862 , primIdTree(MaxPrimId)
1863 , primIdAcc(primIdTree)
1879 if (mPrimCount == MaxPrimId || primIdTree.leafCount() > 1000) {
1884 return mPrimCount++;
1889 enum { MaxPrimId = 100 };
1891 unsigned char mPrimCount;
1895 template<
typename TreeType,
typename MeshDataAdapter,
typename Interrupter = util::NullInterrupter>
1901 typedef tbb::enumerable_thread_specific<typename VoxelizationDataType::Ptr>
DataTable;
1905 Interrupter* interrupter = NULL)
1906 : mDataTable(&dataTable)
1908 , mInterrupter(interrupter)
1915 if (!dataPtr) dataPtr.reset(
new VoxelizationDataType());
1919 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1922 tbb::task::self().cancel_group_execution();
1926 const size_t numVerts = mMesh->vertexCount(n);
1929 if (numVerts == 3 || numVerts == 4) {
1931 prim.index =
Int32(n);
1933 mMesh->getIndexSpacePoint(n, 0, prim.a);
1934 mMesh->getIndexSpacePoint(n, 1, prim.b);
1935 mMesh->getIndexSpacePoint(n, 2, prim.c);
1937 evalTriangle(prim, *dataPtr);
1939 if (numVerts == 4) {
1940 mMesh->getIndexSpacePoint(n, 3, prim.b);
1941 evalTriangle(prim, *dataPtr);
1949 bool wasInterrupted()
const {
return mInterrupter && mInterrupter->wasInterrupted(); }
1951 struct Triangle {
Vec3d a, b, c;
Int32 index; };
1955 SubTask(
const Triangle& prim, DataTable& dataTable,
size_t polygonCount)
1956 : mLocalDataTable(&dataTable)
1958 , mPolygonCount(polygonCount)
1962 void operator()()
const
1964 const size_t minNumTask = size_t(tbb::task_scheduler_init::default_num_threads() * 10);
1966 if (mPolygonCount > minNumTask) {
1968 typename VoxelizationDataType::Ptr& dataPtr = mLocalDataTable->local();
1969 if (!dataPtr) dataPtr.reset(
new VoxelizationDataType());
1971 voxelizeTriangle(mPrim, *dataPtr);
1974 spawnTasks(mPrim, *mLocalDataTable, mPolygonCount);
1978 DataTable *
const mLocalDataTable;
1979 const Triangle mPrim;
1980 const size_t mPolygonCount;
1984 void evalTriangle(
const Triangle& prim, VoxelizationDataType& data)
const
1986 const size_t minNumTask = size_t(tbb::task_scheduler_init::default_num_threads() * 10);
1988 if (mMesh->polygonCount() > minNumTask) {
1990 voxelizeTriangle(prim, data);
1994 spawnTasks(prim, *mDataTable, mMesh->polygonCount());
1998 static void spawnTasks(
const Triangle& mainPrim, DataTable& dataTable,
size_t primCount)
2000 const size_t newPrimCount = primCount * 4;
2002 tbb::task_group tasks;
2004 const Vec3d ac = (mainPrim.a + mainPrim.c) * 0.5;
2005 const Vec3d bc = (mainPrim.b + mainPrim.c) * 0.5;
2006 const Vec3d ab = (mainPrim.a + mainPrim.b) * 0.5;
2009 prim.index = mainPrim.index;
2011 prim.a = mainPrim.a;
2014 tasks.run(SubTask(prim, dataTable, newPrimCount));
2019 tasks.run(SubTask(prim, dataTable, newPrimCount));
2022 prim.b = mainPrim.b;
2024 tasks.run(SubTask(prim, dataTable, newPrimCount));
2028 prim.c = mainPrim.c;
2029 tasks.run(SubTask(prim, dataTable, newPrimCount));
2034 static void voxelizeTriangle(
const Triangle& prim, VoxelizationDataType& data)
2036 std::deque<Coord> coordList;
2039 ijk = Coord::floor(prim.a);
2040 coordList.push_back(ijk);
2042 computeDistance(ijk, prim, data);
2044 unsigned char primId = data.getNewPrimId();
2045 data.primIdAcc.setValueOnly(ijk, primId);
2047 while (!coordList.empty()) {
2048 ijk = coordList.back();
2049 coordList.pop_back();
2051 for (
Int32 i = 0; i < 26; ++i) {
2053 if (primId != data.primIdAcc.getValue(nijk)) {
2054 data.primIdAcc.setValueOnly(nijk, primId);
2055 if(computeDistance(nijk, prim, data)) coordList.push_back(nijk);
2061 static bool computeDistance(
const Coord& ijk,
const Triangle& prim, VoxelizationDataType& data)
2063 Vec3d uvw, voxelCenter(ijk[0], ijk[1], ijk[2]);
2065 typedef typename TreeType::ValueType ValueType;
2067 const ValueType dist = ValueType((voxelCenter -
2070 const ValueType oldDist = data.distAcc.getValue(ijk);
2072 if (dist < oldDist) {
2073 data.distAcc.setValue(ijk, dist);
2074 data.indexAcc.setValue(ijk, prim.index);
2078 data.indexAcc.setValueOnly(ijk,
std::min(prim.index, data.indexAcc.getValue(ijk)));
2081 return !(dist > 0.75);
2084 DataTable *
const mDataTable;
2086 Interrupter *
const mInterrupter;
2093 template<
typename TreeType>
2099 typedef typename TreeType::template ValueConverter<bool>::Type
BoolTreeType;
2103 std::vector<BoolLeafNodeType*>& lhsNodes)
2104 : mRhsTree(&rhsTree), mLhsNodes(lhsNodes.empty() ? NULL : &lhsNodes[0])
2112 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2114 BoolLeafNodeType* lhsNode = mLhsNodes[n];
2115 const LeafNodeType* rhsNode = acc.
probeConstLeaf(lhsNode->origin());
2117 if (rhsNode) lhsNode->topologyDifference(*rhsNode,
false);
2122 TreeType
const *
const mRhsTree;
2123 BoolLeafNodeType **
const mLhsNodes;
2127 template<
typename LeafNodeTypeA,
typename LeafNodeTypeB>
2131 : mNodesA(nodesA.empty() ? NULL : &nodesA[0])
2132 , mNodesB(nodesB.empty() ? NULL : &nodesB[0])
2137 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2138 mNodesA[n]->topologyUnion(*mNodesB[n]);
2143 LeafNodeTypeA **
const mNodesA;
2144 LeafNodeTypeB **
const mNodesB;
2148 template<
typename TreeType>
2153 typedef typename TreeType::template ValueConverter<bool>::Type
BoolTreeType;
2158 , mNodes(nodes.empty() ? NULL : &nodes[0])
2159 , mLocalMaskTree(false)
2160 , mMaskTree(&maskTree)
2166 , mNodes(rhs.mNodes)
2167 , mLocalMaskTree(false)
2168 , mMaskTree(&mLocalMaskTree)
2174 typedef typename LeafNodeType::ValueOnCIter Iterator;
2179 Coord ijk, nijk, localCorod;
2182 for (
size_t n = range.begin(); n != range.end(); ++n) {
2184 LeafNodeType& node = *mNodes[n];
2186 CoordBBox bbox = node.getNodeBoundingBox();
2189 BoolLeafNodeType& maskNode = *maskAcc.
touchLeaf(node.origin());
2191 for (Iterator it = node.cbeginValueOn(); it; ++it) {
2192 ijk = it.getCoord();
2195 localCorod = LeafNodeType::offsetToLocalCoord(pos);
2197 if (localCorod[2] <
int(LeafNodeType::DIM - 1)) {
2199 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2201 nijk = ijk.offsetBy(0, 0, 1);
2205 if (localCorod[2] > 0) {
2207 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2209 nijk = ijk.offsetBy(0, 0, -1);
2213 if (localCorod[1] <
int(LeafNodeType::DIM - 1)) {
2214 npos = pos + LeafNodeType::DIM;
2215 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2217 nijk = ijk.offsetBy(0, 1, 0);
2221 if (localCorod[1] > 0) {
2222 npos = pos - LeafNodeType::DIM;
2223 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2225 nijk = ijk.offsetBy(0, -1, 0);
2229 if (localCorod[0] <
int(LeafNodeType::DIM - 1)) {
2230 npos = pos + LeafNodeType::DIM * LeafNodeType::DIM;
2231 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2233 nijk = ijk.offsetBy(1, 0, 0);
2237 if (localCorod[0] > 0) {
2238 npos = pos - LeafNodeType::DIM * LeafNodeType::DIM;
2239 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2241 nijk = ijk.offsetBy(-1, 0, 0);
2251 TreeType
const *
const mTree;
2252 LeafNodeType **
const mNodes;
2254 BoolTreeType mLocalMaskTree;
2255 BoolTreeType *
const mMaskTree;
2260 template<
typename TreeType,
typename MeshDataAdapter>
2267 typedef typename TreeType::template ValueConverter<bool>::Type
BoolTreeType;
2271 std::vector<BoolLeafNodeType*>& maskNodes,
2272 BoolTreeType& maskTree,
2274 Int32TreeType& indexTree,
2276 ValueType exteriorBandWidth,
2277 ValueType interiorBandWidth,
2278 ValueType voxelSize)
2279 : mMaskNodes(maskNodes.empty() ? NULL : &maskNodes[0])
2280 , mMaskTree(&maskTree)
2281 , mDistTree(&distTree)
2282 , mIndexTree(&indexTree)
2284 , mNewMaskTree(false)
2286 , mUpdatedDistNodes()
2288 , mUpdatedIndexNodes()
2289 , mExteriorBandWidth(exteriorBandWidth)
2290 , mInteriorBandWidth(interiorBandWidth)
2291 , mVoxelSize(voxelSize)
2296 : mMaskNodes(rhs.mMaskNodes)
2297 , mMaskTree(rhs.mMaskTree)
2298 , mDistTree(rhs.mDistTree)
2299 , mIndexTree(rhs.mIndexTree)
2301 , mNewMaskTree(false)
2303 , mUpdatedDistNodes()
2305 , mUpdatedIndexNodes()
2306 , mExteriorBandWidth(rhs.mExteriorBandWidth)
2307 , mInteriorBandWidth(rhs.mInteriorBandWidth)
2308 , mVoxelSize(rhs.mVoxelSize)
2313 mDistNodes.insert(mDistNodes.end(), rhs.mDistNodes.begin(), rhs.mDistNodes.end());
2314 mIndexNodes.insert(mIndexNodes.end(), rhs.mIndexNodes.begin(), rhs.mIndexNodes.end());
2316 mUpdatedDistNodes.insert(mUpdatedDistNodes.end(), rhs.mUpdatedDistNodes.begin(), rhs.mUpdatedDistNodes.end());
2317 mUpdatedIndexNodes.insert(mUpdatedIndexNodes.end(), rhs.mUpdatedIndexNodes.begin(), rhs.mUpdatedIndexNodes.end());
2319 mNewMaskTree.merge(rhs.mNewMaskTree);
2329 std::vector<Int32> primitives;
2330 primitives.reserve(26);
2332 LeafNodeType * newDistNodePt = NULL;
2333 Int32LeafNodeType * newIndexNodePt = NULL;
2335 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2337 BoolLeafNodeType& maskNode = *mMaskNodes[n];
2338 if (maskNode.isEmpty())
continue;
2340 Coord ijk = maskNode.origin();
2342 bool usingNewNodes =
false;
2344 LeafNodeType * distNodePt = distAcc.
probeLeaf(ijk);
2345 Int32LeafNodeType * indexNodePt = indexAcc.
probeLeaf(ijk);
2347 assert(!distNodePt == !indexNodePt);
2349 if (!distNodePt && !indexNodePt) {
2351 const ValueType backgroundDist = distAcc.
getValue(ijk);
2353 if (!newDistNodePt && !newIndexNodePt) {
2354 newDistNodePt =
new LeafNodeType(ijk, backgroundDist);
2355 newIndexNodePt =
new Int32LeafNodeType(ijk, indexAcc.
getValue(ijk));
2358 if ((backgroundDist < ValueType(0.0)) != (newDistNodePt->getValue(0) < ValueType(0.0))) {
2359 newDistNodePt->buffer().fill(backgroundDist);
2362 newDistNodePt->setOrigin(ijk);
2363 newIndexNodePt->setOrigin(ijk);
2366 distNodePt = newDistNodePt;
2367 indexNodePt = newIndexNodePt;
2369 usingNewNodes =
true;
2372 bool updatedValues =
false;
2374 for (
typename BoolLeafNodeType::ValueOnIter it = maskNode.beginValueOn(); it; ++it) {
2376 ijk = it.getCoord();
2377 const Index pos = it.pos();
2379 Int32 closestPrimIdx = 0;
2380 const ValueType distance =
2381 computeDistance(ijk, distAcc, indexAcc, primitives, closestPrimIdx);
2383 const bool inside = distNodePt->getValue(pos) < ValueType(0.0);
2385 if (!inside && distance < mExteriorBandWidth) {
2386 distNodePt->setValueOnly(pos, distance);
2387 indexNodePt->setValueOn(pos, closestPrimIdx);
2388 }
else if (inside && distance < mInteriorBandWidth) {
2389 distNodePt->setValueOnly(pos, -distance);
2390 indexNodePt->setValueOn(pos, closestPrimIdx);
2395 for (
Int32 i = 0; i < 6; ++i) {
2396 newMaskAcc.
setValueOn(ijk + util::COORD_OFFSETS[i]);
2399 updatedValues =
true;
2403 if (updatedValues && usingNewNodes) {
2405 distNodePt->topologyUnion(*indexNodePt);
2407 mDistNodes.push_back(distNodePt);
2408 mIndexNodes.push_back(indexNodePt);
2410 newDistNodePt = NULL;
2411 newIndexNodePt = NULL;
2413 }
else if (updatedValues) {
2415 mUpdatedDistNodes.push_back(distNodePt);
2416 mUpdatedIndexNodes.push_back(indexNodePt);
2436 computeDistance(
const Coord& ijk,
2438 std::vector<Int32>& primitives,
Int32& closestPrimIdx)
const
2443 const Coord ijkMin = ijk.offsetBy(-1);
2444 const Coord ijkMax = ijk.offsetBy(1);
2445 const Coord nodeMin = ijkMin & ~(LeafNodeType::DIM - 1);
2446 const Coord nodeMax = ijkMax & ~(LeafNodeType::DIM - 1);
2451 for (nijk[0] = nodeMin[0]; nijk[0] <= nodeMax[0]; nijk[0] += LeafNodeType::DIM) {
2452 for (nijk[1] = nodeMin[1]; nijk[1] <= nodeMax[1]; nijk[1] += LeafNodeType::DIM) {
2453 for (nijk[2] = nodeMin[2]; nijk[2] <= nodeMax[2]; nijk[2] += LeafNodeType::DIM) {
2455 if (LeafNodeType* distleaf = distAcc.
probeLeaf(nijk)) {
2460 evalLeafNode(bbox, *distleaf, *idxAcc.
probeLeaf(nijk), primitives, minDist);
2466 const ValueType tmpDist = evalPrimitives(ijk, primitives, closestPrimIdx);
2467 return tmpDist > minDist ? tmpDist : minDist + mVoxelSize;
2471 evalLeafNode(
const CoordBBox& bbox, LeafNodeType& distLeaf,
2472 Int32LeafNodeType& idxLeaf, std::vector<Int32>& primitives, ValueType& minNeighbourDist)
const
2475 Index xPos(0), yPos(0), pos(0);
2477 for (
int x = bbox.min()[0]; x <= bbox.max()[0]; ++x) {
2478 xPos = (x & (LeafNodeType::DIM - 1u)) << (2 * LeafNodeType::LOG2DIM);
2479 for (
int y = bbox.min()[1]; y <= bbox.max()[1]; ++y) {
2480 yPos = xPos + ((y & (LeafNodeType::DIM - 1u)) << LeafNodeType::LOG2DIM);
2481 for (
int z = bbox.min()[2]; z <= bbox.max()[2]; ++z) {
2482 pos = yPos + (z & (LeafNodeType::DIM - 1u));
2483 if (distLeaf.probeValue(pos, tmpDist)) {
2484 primitives.push_back(idxLeaf.getValue(pos));
2485 minNeighbourDist =
std::min(std::abs(tmpDist), minNeighbourDist);
2493 evalPrimitives(
const Coord& ijk, std::vector<Int32>& primitives,
Int32& closestPrimIdx)
const
2495 std::sort(primitives.begin(), primitives.end());
2497 Int32 lastPrim = -1;
2498 Vec3d a, b, c, uvw, voxelCenter(ijk[0], ijk[1], ijk[2]);
2500 for (
size_t n = 0, N = primitives.size(); n < N; ++n) {
2502 if (primitives[n] == lastPrim)
continue;
2504 lastPrim = primitives[n];
2506 const size_t polygon = size_t(lastPrim);
2508 mMesh->getIndexSpacePoint(polygon, 0, a);
2509 mMesh->getIndexSpacePoint(polygon, 1, b);
2510 mMesh->getIndexSpacePoint(polygon, 2, c);
2512 primDist = (voxelCenter -
2516 if (4 == mMesh->vertexCount(polygon)) {
2518 mMesh->getIndexSpacePoint(polygon, 3, b);
2521 a, b, c, voxelCenter, uvw)).lengthSqr();
2523 if (tmpDist < primDist) primDist = tmpDist;
2526 if (primDist < dist) {
2528 closestPrimIdx = lastPrim;
2532 return ValueType(std::sqrt(dist)) * mVoxelSize;
2539 BoolLeafNodeType **
const mMaskNodes;
2540 BoolTreeType *
const mMaskTree;
2541 TreeType *
const mDistTree;
2542 Int32TreeType *
const mIndexTree;
2546 BoolTreeType mNewMaskTree;
2548 std::vector<LeafNodeType*> mDistNodes, mUpdatedDistNodes;
2549 std::vector<Int32LeafNodeType*> mIndexNodes, mUpdatedIndexNodes;
2551 const ValueType mExteriorBandWidth, mInteriorBandWidth, mVoxelSize;
2555 template<
typename TreeType,
typename Int32TreeType,
typename BoolTreeType,
typename MeshDataAdapter>
2559 Int32TreeType& indexTree,
2560 BoolTreeType& maskTree,
2561 std::vector<typename BoolTreeType::LeafNodeType*>& maskNodes,
2563 typename TreeType::ValueType exteriorBandWidth,
2564 typename TreeType::ValueType interiorBandWidth,
2565 typename TreeType::ValueType voxelSize)
2567 typedef typename TreeType::LeafNodeType LeafNodeType;
2568 typedef typename Int32TreeType::LeafNodeType Int32LeafNodeType;
2571 expandOp(maskNodes, maskTree, distTree, indexTree,
2572 mesh, exteriorBandWidth, interiorBandWidth, voxelSize);
2574 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, maskNodes.size()), expandOp);
2578 typedef typename std::vector<LeafNodeType*> LeafNodePtVec;
2580 for (
typename LeafNodePtVec::iterator it = nodes.begin(), end = nodes.end(); it != end; ++it) {
2587 typedef typename std::vector<Int32LeafNodeType*> LeafNodePtVec;
2589 for (
typename LeafNodePtVec::iterator it = nodes.begin(), end = nodes.end(); it != end; ++it) {
2594 tbb::parallel_for(tbb::blocked_range<size_t>(0, expandOp.
updatedIndexNodes().size()),
2606 template<
typename TreeType>
2613 ValueType voxelSize,
bool unsignedDist)
2615 , mVoxelSize(voxelSize)
2616 , mUnsigned(unsignedDist)
2622 typename LeafNodeType::ValueOnIter iter;
2624 const bool udf = mUnsigned;
2625 const ValueType w[2] = { -mVoxelSize, mVoxelSize };
2627 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2629 for (iter = mNodes[n]->beginValueOn(); iter; ++iter) {
2630 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
2631 val = w[udf || (val < ValueType(0.0))] * std::sqrt(std::abs(val));
2637 LeafNodeType * *
const mNodes;
2638 const ValueType mVoxelSize;
2639 const bool mUnsigned;
2644 template<
typename TreeType>
2651 ValueType exBandWidth, ValueType inBandWidth)
2652 : mNodes(nodes.empty() ? NULL : &nodes[0])
2653 , mExBandWidth(exBandWidth)
2654 , mInBandWidth(inBandWidth)
2660 typename LeafNodeType::ValueOnIter iter;
2661 const ValueType exVal = mExBandWidth;
2662 const ValueType inVal = -mInBandWidth;
2664 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2666 for (iter = mNodes[n]->beginValueOn(); iter; ++iter) {
2668 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
2670 const bool inside = val < ValueType(0.0);
2672 if (inside && !(val > inVal)) {
2675 }
else if (!inside && !(val < exVal)) {
2684 LeafNodeType * *
const mNodes;
2685 const ValueType mExBandWidth, mInBandWidth;
2689 template<
typename TreeType>
2696 : mNodes(nodes.empty() ? NULL : &nodes[0]), mOffset(offset)
2702 const ValueType offset = mOffset;
2704 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2706 typename LeafNodeType::ValueOnIter iter = mNodes[n]->beginValueOn();
2708 for (; iter; ++iter) {
2709 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
2716 LeafNodeType * *
const mNodes;
2717 const ValueType mOffset;
2721 template<
typename TreeType>
2727 Renormalize(
const TreeType& tree,
const std::vector<LeafNodeType*>& nodes, ValueType* buffer, ValueType voxelSize)
2729 , mNodes(nodes.empty() ? NULL : &nodes[0])
2731 , mVoxelSize(voxelSize)
2744 const ValueType dx = mVoxelSize, invDx = ValueType(1.0) / mVoxelSize;
2746 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2748 ValueType* bufferData = &mBuffer[n * LeafNodeType::SIZE];
2750 typename LeafNodeType::ValueOnCIter iter = mNodes[n]->cbeginValueOn();
2751 for (; iter; ++iter) {
2753 const ValueType phi0 = *iter;
2755 ijk = iter.getCoord();
2757 up[0] = acc.
getValue(ijk.offsetBy(1, 0, 0)) - phi0;
2758 up[1] = acc.
getValue(ijk.offsetBy(0, 1, 0)) - phi0;
2759 up[2] = acc.
getValue(ijk.offsetBy(0, 0, 1)) - phi0;
2761 down[0] = phi0 - acc.
getValue(ijk.offsetBy(-1, 0, 0));
2762 down[1] = phi0 - acc.
getValue(ijk.offsetBy(0, -1, 0));
2763 down[2] = phi0 - acc.
getValue(ijk.offsetBy(0, 0, -1));
2767 const ValueType diff =
math::Sqrt(normSqGradPhi) * invDx - ValueType(1.0);
2770 bufferData[iter.pos()] = phi0 - dx * S * diff;
2776 TreeType
const *
const mTree;
2777 LeafNodeType
const *
const *
const mNodes;
2778 ValueType *
const mBuffer;
2780 const ValueType mVoxelSize;
2784 template<
typename TreeType>
2790 MinCombine(std::vector<LeafNodeType*>& nodes,
const ValueType* buffer)
2791 : mNodes(nodes.empty() ? NULL : &nodes[0]), mBuffer(buffer)
2797 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2799 const ValueType* bufferData = &mBuffer[n * LeafNodeType::SIZE];
2801 typename LeafNodeType::ValueOnIter iter = mNodes[n]->beginValueOn();
2803 for (; iter; ++iter) {
2804 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
2805 val =
std::min(val, bufferData[iter.pos()]);
2811 LeafNodeType * *
const mNodes;
2812 ValueType
const *
const mBuffer;
2824 template <
typename FloatTreeT>
2830 ConnectivityTable nodeConnectivity(tree);
2832 std::vector<size_t> zStartNodes, yStartNodes, xStartNodes;
2834 for (
size_t n = 0; n < nodeConnectivity.size(); ++n) {
2835 if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevX()[n]) {
2836 xStartNodes.push_back(n);
2839 if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevY()[n]) {
2840 yStartNodes.push_back(n);
2843 if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevZ()[n]) {
2844 zStartNodes.push_back(n);
2850 tbb::parallel_for(tbb::blocked_range<size_t>(0, zStartNodes.size()),
2853 tbb::parallel_for(tbb::blocked_range<size_t>(0, yStartNodes.size()),
2856 tbb::parallel_for(tbb::blocked_range<size_t>(0, xStartNodes.size()),
2859 const size_t numLeafNodes = nodeConnectivity.size();
2860 const size_t numVoxels = numLeafNodes * FloatTreeT::LeafNodeType::SIZE;
2862 boost::scoped_array<bool> changedNodeMaskA(
new bool[numLeafNodes]);
2863 boost::scoped_array<bool> changedNodeMaskB(
new bool[numLeafNodes]);
2864 boost::scoped_array<bool> changedVoxelMask(
new bool[numVoxels]);
2866 memset(changedNodeMaskA.get(), 1,
sizeof(bool) * numLeafNodes);
2869 const tbb::blocked_range<size_t> nodeRange(0, numLeafNodes);
2871 bool nodesUpdated =
false;
2874 nodeConnectivity.nodes(), changedNodeMaskA.get()));
2877 changedNodeMaskA.get(), changedNodeMaskB.get(), changedVoxelMask.get()));
2879 changedNodeMaskA.swap(changedNodeMaskB);
2881 nodesUpdated =
false;
2882 for (
size_t n = 0; n < numLeafNodes; ++n) {
2883 nodesUpdated |= changedNodeMaskA[n];
2884 if (nodesUpdated)
break;
2889 nodeConnectivity.nodes(), changedNodeMaskA.get(), changedVoxelMask.get()));
2891 }
while (nodesUpdated);
2899 template <
typename Gr
idType,
typename MeshDataAdapter,
typename Interrupter>
2900 inline typename GridType::Ptr
2902 Interrupter& interrupter,
2905 float exteriorBandWidth,
2906 float interiorBandWidth,
2910 typedef typename GridType::Ptr GridTypePtr;
2911 typedef typename GridType::TreeType TreeType;
2912 typedef typename TreeType::LeafNodeType LeafNodeType;
2913 typedef typename GridType::ValueType ValueType;
2916 typedef typename Int32GridType::TreeType Int32TreeType;
2918 typedef typename TreeType::template ValueConverter<bool>::Type BoolTreeType;
2925 distGrid->setTransform(transform.
copy());
2927 ValueType exteriorWidth = ValueType(exteriorBandWidth);
2928 ValueType interiorWidth = ValueType(interiorBandWidth);
2932 if (!boost::math::isfinite(exteriorWidth) || boost::math::isnan(interiorWidth)) {
2933 std::stringstream msg;
2934 msg <<
"Illegal narrow band width: exterior = " << exteriorWidth
2935 <<
", interior = " << interiorWidth;
2940 const ValueType voxelSize = ValueType(transform.
voxelSize()[0]);
2942 if (!boost::math::isfinite(voxelSize) ||
math::isZero(voxelSize)) {
2943 std::stringstream msg;
2944 msg <<
"Illegal transform, voxel size = " << voxelSize;
2950 exteriorWidth *= voxelSize;
2954 interiorWidth *= voxelSize;
2962 Int32GridType* indexGrid = NULL;
2964 typename Int32GridType::Ptr temporaryIndexGrid;
2966 if (polygonIndexGrid) {
2967 indexGrid = polygonIndexGrid;
2970 indexGrid = temporaryIndexGrid.get();
2973 indexGrid->newTree();
2974 indexGrid->setTransform(transform.
copy());
2976 if (computeSignedDistanceField) {
2980 interiorWidth = ValueType(0.0);
2983 TreeType& distTree = distGrid->tree();
2984 Int32TreeType& indexTree = indexGrid->tree();
2993 typedef tbb::enumerable_thread_specific<typename VoxelizationDataType::Ptr> DataTable;
2998 const tbb::blocked_range<size_t> polygonRange(0, mesh.polygonCount());
3000 tbb::parallel_for(polygonRange, Voxelizer(data, mesh, &interrupter));
3002 for (
typename DataTable::iterator i = data.begin(); i != data.end(); ++i) {
3003 VoxelizationDataType& dataItem = **i;
3010 if (interrupter.wasInterrupted(30))
return distGrid;
3017 if (computeSignedDistanceField) {
3022 std::vector<LeafNodeType*> nodes;
3023 nodes.reserve(distTree.leafCount());
3024 distTree.getNodes(nodes);
3026 const tbb::blocked_range<size_t> nodeRange(0, nodes.size());
3030 tbb::parallel_for(nodeRange, SignOp(nodes, distTree, indexTree, mesh));
3032 if (interrupter.wasInterrupted(45))
return distGrid;
3035 if (removeIntersectingVoxels) {
3037 tbb::parallel_for(nodeRange,
3040 tbb::parallel_for(nodeRange,
3048 if (interrupter.wasInterrupted(50))
return distGrid;
3050 if (distTree.activeVoxelCount() == 0) {
3051 distGrid.reset((
new GridType(ValueType(0.0))));
3057 std::vector<LeafNodeType*> nodes;
3058 nodes.reserve(distTree.leafCount());
3059 distTree.getNodes(nodes);
3061 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3066 if (computeSignedDistanceField) {
3067 distTree.root().setBackground(exteriorWidth,
false);
3073 if (interrupter.wasInterrupted(54))
return distGrid;
3080 const ValueType minBandWidth = voxelSize * ValueType(2.0);
3082 if (interiorWidth > minBandWidth || exteriorWidth > minBandWidth) {
3085 BoolTreeType maskTree(
false);
3088 std::vector<LeafNodeType*> nodes;
3089 nodes.reserve(distTree.leafCount());
3090 distTree.getNodes(nodes);
3093 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), op);
3099 float progress = 54.0f, step = 0.0f;
3101 2.0 * std::ceil((
std::max(interiorWidth, exteriorWidth) - minBandWidth) / voxelSize);
3103 if (estimated <
double(maxIterations)) {
3104 maxIterations = unsigned(estimated);
3105 step = 40.0f / float(maxIterations);
3108 std::vector<typename BoolTreeType::LeafNodeType*> maskNodes;
3114 if (interrupter.wasInterrupted(
int(progress)))
return distGrid;
3116 const size_t maskNodeCount = maskTree.leafCount();
3117 if (maskNodeCount == 0)
break;
3120 maskNodes.reserve(maskNodeCount);
3121 maskTree.getNodes(maskNodes);
3123 const tbb::blocked_range<size_t> range(0, maskNodes.size());
3125 tbb::parallel_for(range,
3129 mesh, exteriorWidth, interiorWidth, voxelSize);
3131 if ((++count) >= maxIterations)
break;
3136 if (interrupter.wasInterrupted(94))
return distGrid;
3138 if (!polygonIndexGrid) indexGrid->clear();
3146 if (computeSignedDistanceField && renormalizeValues) {
3148 std::vector<LeafNodeType*> nodes;
3149 nodes.reserve(distTree.leafCount());
3150 distTree.getNodes(nodes);
3152 boost::scoped_array<ValueType> buffer(
new ValueType[LeafNodeType::SIZE * nodes.size()]);
3154 const ValueType offset = ValueType(0.8 * voxelSize);
3156 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3159 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3162 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3165 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3169 if (interrupter.wasInterrupted(99))
return distGrid;
3176 if (trimNarrowBand &&
std::min(interiorWidth, exteriorWidth) < voxelSize * ValueType(4.0)) {
3178 std::vector<LeafNodeType*> nodes;
3179 nodes.reserve(distTree.leafCount());
3180 distTree.getNodes(nodes);
3182 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3185 tools::pruneLevelSet(distTree, exteriorWidth, computeSignedDistanceField ? -interiorWidth : -exteriorWidth);
3192 template <
typename Gr
idType,
typename MeshDataAdapter>
3193 inline typename GridType::Ptr
3197 float exteriorBandWidth,
3198 float interiorBandWidth,
3203 return meshToVolume<GridType>(nullInterrupter, mesh, transform,
3204 exteriorBandWidth, interiorBandWidth, flags, polygonIndexGrid);
3212 template<
typename Gr
idType>
3213 inline typename boost::enable_if<boost::is_floating_point<typename GridType::ValueType>,
3214 typename GridType::Ptr>::type
3216 const openvdb::math::Transform& xform,
3217 const std::vector<Vec3s>& points,
3218 const std::vector<Vec3I>& triangles,
3219 const std::vector<Vec4I>& quads,
3222 bool unsignedDistanceField =
false)
3224 if (points.empty()) {
3225 return typename GridType::Ptr(
new GridType(
typename GridType::ValueType(exBandWidth)));
3228 const size_t numPoints = points.size();
3229 boost::scoped_array<Vec3s> indexSpacePoints(
new Vec3s[numPoints]);
3232 tbb::parallel_for(tbb::blocked_range<size_t>(0, numPoints),
3234 &points[0], indexSpacePoints.get(), xform));
3238 if (quads.empty()) {
3241 mesh(indexSpacePoints.get(), numPoints, &triangles[0], triangles.size());
3243 return meshToVolume<GridType>(mesh, xform, exBandWidth, inBandWidth, conversionFlags);
3245 }
else if (triangles.empty()) {
3248 mesh(indexSpacePoints.get(), numPoints, &quads[0], quads.size());
3250 return meshToVolume<GridType>(mesh, xform, exBandWidth, inBandWidth, conversionFlags);
3255 const size_t numPrimitives = triangles.size() + quads.size();
3256 boost::scoped_array<Vec4I> prims(
new Vec4I[numPrimitives]);
3258 for (
size_t n = 0, N = triangles.size(); n < N; ++n) {
3259 const Vec3I& triangle = triangles[n];
3260 Vec4I& prim = prims[n];
3261 prim[0] = triangle[0];
3262 prim[1] = triangle[1];
3263 prim[2] = triangle[2];
3267 const size_t offset = triangles.size();
3268 for (
size_t n = 0, N = quads.size(); n < N; ++n) {
3269 prims[offset + n] = quads[n];
3273 mesh(indexSpacePoints.get(), numPoints, prims.get(), numPrimitives);
3275 return meshToVolume<GridType>(mesh, xform, exBandWidth, inBandWidth, conversionFlags);
3281 template<
typename Gr
idType>
3282 inline typename boost::disable_if<boost::is_floating_point<typename GridType::ValueType>,
3283 typename GridType::Ptr>::type
3286 const std::vector<Vec3s>& ,
3287 const std::vector<Vec3I>& ,
3288 const std::vector<Vec4I>& ,
3294 "mesh to volume conversion is supported only for scalar floating-point grids");
3301 template<
typename Gr
idType>
3302 inline typename GridType::Ptr
3304 const openvdb::math::Transform& xform,
3305 const std::vector<Vec3s>& points,
3306 const std::vector<Vec3I>& triangles,
3309 std::vector<Vec4I> quads(0);
3310 return doMeshConversion<GridType>(xform, points, triangles, quads,
3311 halfWidth, halfWidth);
3315 template<
typename Gr
idType>
3316 inline typename GridType::Ptr
3318 const openvdb::math::Transform& xform,
3319 const std::vector<Vec3s>& points,
3320 const std::vector<Vec4I>& quads,
3323 std::vector<Vec3I> triangles(0);
3324 return doMeshConversion<GridType>(xform, points, triangles, quads,
3325 halfWidth, halfWidth);
3329 template<
typename Gr
idType>
3330 inline typename GridType::Ptr
3332 const openvdb::math::Transform& xform,
3333 const std::vector<Vec3s>& points,
3334 const std::vector<Vec3I>& triangles,
3335 const std::vector<Vec4I>& quads,
3338 return doMeshConversion<GridType>(xform, points, triangles, quads,
3339 halfWidth, halfWidth);
3343 template<
typename Gr
idType>
3344 inline typename GridType::Ptr
3346 const openvdb::math::Transform& xform,
3347 const std::vector<Vec3s>& points,
3348 const std::vector<Vec3I>& triangles,
3349 const std::vector<Vec4I>& quads,
3353 return doMeshConversion<GridType>(xform, points, triangles,
3354 quads, exBandWidth, inBandWidth);
3358 template<
typename Gr
idType>
3359 inline typename GridType::Ptr
3361 const openvdb::math::Transform& xform,
3362 const std::vector<Vec3s>& points,
3363 const std::vector<Vec3I>& triangles,
3364 const std::vector<Vec4I>& quads,
3367 return doMeshConversion<GridType>(xform, points, triangles, quads,
3368 bandWidth, bandWidth,
true);
3376 inline std::ostream&
3379 ostr <<
"{[ " << rhs.
mXPrim <<
", " << rhs.
mXDist <<
"]";
3380 ostr <<
" [ " << rhs.
mYPrim <<
", " << rhs.
mYDist <<
"]";
3381 ostr <<
" [ " << rhs.
mZPrim <<
", " << rhs.
mZDist <<
"]}";
3386 inline MeshToVoxelEdgeData::EdgeData
3401 const std::vector<Vec3s>& pointList,
3402 const std::vector<Vec4I>& polygonList);
3404 void run(
bool threaded =
true);
3407 inline void operator() (
const tbb::blocked_range<size_t> &range);
3415 struct Primitive {
Vec3d a, b, c, d;
Int32 index; };
3417 template<
bool IsQuad>
3418 inline void voxelize(
const Primitive&);
3420 template<
bool IsQuad>
3421 inline bool evalPrimitive(
const Coord&,
const Primitive&);
3423 inline bool rayTriangleIntersection(
const Vec3d& origin,
const Vec3d& dir,
3430 const std::vector<Vec3s>& mPointList;
3431 const std::vector<Vec4I>& mPolygonList;
3434 typedef TreeType::ValueConverter<Int32>::Type IntTreeT;
3435 IntTreeT mLastPrimTree;
3441 MeshToVoxelEdgeData::GenEdgeData::GenEdgeData(
3442 const std::vector<Vec3s>& pointList,
3443 const std::vector<Vec4I>& polygonList)
3446 , mPointList(pointList)
3447 , mPolygonList(polygonList)
3449 , mLastPrimAccessor(mLastPrimTree)
3458 , mPointList(rhs.mPointList)
3459 , mPolygonList(rhs.mPolygonList)
3461 , mLastPrimAccessor(mLastPrimTree)
3470 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPolygonList.size()), *
this);
3472 (*this)(tbb::blocked_range<size_t>(0, mPolygonList.size()));
3481 typedef RootNodeType::NodeChainType NodeChainType;
3482 BOOST_STATIC_ASSERT(boost::mpl::size<NodeChainType>::value > 1);
3483 typedef boost::mpl::at<NodeChainType, boost::mpl::int_<1> >::type InternalNodeType;
3491 for ( ; leafIt; ++leafIt) {
3492 ijk = leafIt->origin();
3498 mAccessor.addLeaf(rhs.mAccessor.
probeLeaf(ijk));
3499 InternalNodeType* node = rhs.mAccessor.
getNode<InternalNodeType>();
3501 rhs.mAccessor.
clear();
3511 if (!lhsLeafPt->isValueOn(offset)) {
3512 lhsLeafPt->setValueOn(offset, rhsValue);
3544 for (
size_t n = range.begin(); n < range.end(); ++n) {
3546 const Vec4I& verts = mPolygonList[n];
3548 prim.index =
Int32(n);
3549 prim.a =
Vec3d(mPointList[verts[0]]);
3550 prim.b =
Vec3d(mPointList[verts[1]]);
3551 prim.c =
Vec3d(mPointList[verts[2]]);
3554 prim.d =
Vec3d(mPointList[verts[3]]);
3555 voxelize<true>(prim);
3557 voxelize<false>(prim);
3563 template<
bool IsQuad>
3565 MeshToVoxelEdgeData::GenEdgeData::voxelize(
const Primitive& prim)
3567 std::deque<Coord> coordList;
3570 ijk = Coord::floor(prim.a);
3571 coordList.push_back(ijk);
3573 evalPrimitive<IsQuad>(ijk, prim);
3575 while (!coordList.empty()) {
3577 ijk = coordList.back();
3578 coordList.pop_back();
3580 for (
Int32 i = 0; i < 26; ++i) {
3581 nijk = ijk + util::COORD_OFFSETS[i];
3583 if (prim.index != mLastPrimAccessor.getValue(nijk)) {
3584 mLastPrimAccessor.setValue(nijk, prim.index);
3585 if(evalPrimitive<IsQuad>(nijk, prim)) coordList.push_back(nijk);
3592 template<
bool IsQuad>
3594 MeshToVoxelEdgeData::GenEdgeData::evalPrimitive(
const Coord& ijk,
const Primitive& prim)
3596 Vec3d uvw, org(ijk[0], ijk[1], ijk[2]);
3597 bool intersecting =
false;
3601 mAccessor.probeValue(ijk, edgeData);
3604 double dist = (org -
3607 if (rayTriangleIntersection(org,
Vec3d(1.0, 0.0, 0.0), prim.a, prim.c, prim.b, t)) {
3608 if (t < edgeData.mXDist) {
3609 edgeData.mXDist = float(t);
3610 edgeData.mXPrim = prim.index;
3611 intersecting =
true;
3615 if (rayTriangleIntersection(org,
Vec3d(0.0, 1.0, 0.0), prim.a, prim.c, prim.b, t)) {
3616 if (t < edgeData.mYDist) {
3617 edgeData.mYDist = float(t);
3618 edgeData.mYPrim = prim.index;
3619 intersecting =
true;
3623 if (rayTriangleIntersection(org,
Vec3d(0.0, 0.0, 1.0), prim.a, prim.c, prim.b, t)) {
3624 if (t < edgeData.mZDist) {
3625 edgeData.mZDist = float(t);
3626 edgeData.mZPrim = prim.index;
3627 intersecting =
true;
3633 double secondDist = (org -
3636 if (secondDist < dist) dist = secondDist;
3638 if (rayTriangleIntersection(org,
Vec3d(1.0, 0.0, 0.0), prim.a, prim.d, prim.c, t)) {
3639 if (t < edgeData.mXDist) {
3640 edgeData.mXDist = float(t);
3641 edgeData.mXPrim = prim.index;
3642 intersecting =
true;
3646 if (rayTriangleIntersection(org,
Vec3d(0.0, 1.0, 0.0), prim.a, prim.d, prim.c, t)) {
3647 if (t < edgeData.mYDist) {
3648 edgeData.mYDist = float(t);
3649 edgeData.mYPrim = prim.index;
3650 intersecting =
true;
3654 if (rayTriangleIntersection(org,
Vec3d(0.0, 0.0, 1.0), prim.a, prim.d, prim.c, t)) {
3655 if (t < edgeData.mZDist) {
3656 edgeData.mZDist = float(t);
3657 edgeData.mZPrim = prim.index;
3658 intersecting =
true;
3663 if (intersecting) mAccessor.setValue(ijk, edgeData);
3665 return (dist < 0.86602540378443861);
3670 MeshToVoxelEdgeData::GenEdgeData::rayTriangleIntersection(
3681 double divisor = s1.
dot(e1);
3682 if (!(std::abs(divisor) > 0.0))
return false;
3686 double inv_divisor = 1.0 / divisor;
3687 Vec3d d = origin - a;
3688 double b1 = d.
dot(s1) * inv_divisor;
3690 if (b1 < 0.0 || b1 > 1.0)
return false;
3693 double b2 = dir.
dot(s2) * inv_divisor;
3695 if (b2 < 0.0 || (b1 + b2) > 1.0)
return false;
3699 t = e2.dot(s2) * inv_divisor;
3700 return (t < 0.0) ?
false :
true;
3716 const std::vector<Vec3s>& pointList,
3717 const std::vector<Vec4I>& polygonList)
3731 std::vector<Vec3d>& points,
3732 std::vector<Index32>& primitives)
3742 point[0] = double(coord[0]) + data.
mXDist;
3743 point[1] = double(coord[1]);
3744 point[2] = double(coord[2]);
3746 points.push_back(point);
3747 primitives.push_back(data.
mXPrim);
3751 point[0] = double(coord[0]);
3752 point[1] = double(coord[1]) + data.
mYDist;
3753 point[2] = double(coord[2]);
3755 points.push_back(point);
3756 primitives.push_back(data.
mYPrim);
3760 point[0] = double(coord[0]);
3761 point[1] = double(coord[1]);
3762 point[2] = double(coord[2]) + data.
mZDist;
3764 points.push_back(point);
3765 primitives.push_back(data.
mZPrim);
3775 point[0] = double(coord[0]);
3776 point[1] = double(coord[1]) + data.
mYDist;
3777 point[2] = double(coord[2]);
3779 points.push_back(point);
3780 primitives.push_back(data.
mYPrim);
3784 point[0] = double(coord[0]);
3785 point[1] = double(coord[1]);
3786 point[2] = double(coord[2]) + data.
mZDist;
3788 points.push_back(point);
3789 primitives.push_back(data.
mZPrim);
3797 point[0] = double(coord[0]);
3798 point[1] = double(coord[1]) + data.
mYDist;
3799 point[2] = double(coord[2]);
3801 points.push_back(point);
3802 primitives.push_back(data.
mYPrim);
3811 point[0] = double(coord[0]) + data.
mXDist;
3812 point[1] = double(coord[1]);
3813 point[2] = double(coord[2]);
3815 points.push_back(point);
3816 primitives.push_back(data.
mXPrim);
3820 point[0] = double(coord[0]);
3821 point[1] = double(coord[1]) + data.
mYDist;
3822 point[2] = double(coord[2]);
3824 points.push_back(point);
3825 primitives.push_back(data.
mYPrim);
3835 point[0] = double(coord[0]) + data.
mXDist;
3836 point[1] = double(coord[1]);
3837 point[2] = double(coord[2]);
3839 points.push_back(point);
3840 primitives.push_back(data.
mXPrim);
3849 point[0] = double(coord[0]) + data.
mXDist;
3850 point[1] = double(coord[1]);
3851 point[2] = double(coord[2]);
3853 points.push_back(point);
3854 primitives.push_back(data.
mXPrim);
3858 point[0] = double(coord[0]);
3859 point[1] = double(coord[1]);
3860 point[2] = double(coord[2]) + data.
mZDist;
3862 points.push_back(point);
3863 primitives.push_back(data.
mZPrim);
3872 point[0] = double(coord[0]);
3873 point[1] = double(coord[1]);
3874 point[2] = double(coord[2]) + data.
mZDist;
3876 points.push_back(point);
3877 primitives.push_back(data.
mZPrim);
3883 template<
typename Gr
idType,
typename VecType>
3884 inline typename GridType::Ptr
3886 const openvdb::math::Transform& xform,
3887 typename VecType::ValueType halfWidth)
3893 points[0] =
Vec3s(pmin[0], pmin[1], pmin[2]);
3894 points[1] =
Vec3s(pmin[0], pmin[1], pmax[2]);
3895 points[2] =
Vec3s(pmax[0], pmin[1], pmax[2]);
3896 points[3] =
Vec3s(pmax[0], pmin[1], pmin[2]);
3897 points[4] =
Vec3s(pmin[0], pmax[1], pmin[2]);
3898 points[5] =
Vec3s(pmin[0], pmax[1], pmax[2]);
3899 points[6] =
Vec3s(pmax[0], pmax[1], pmax[2]);
3900 points[7] =
Vec3s(pmax[0], pmax[1], pmin[2]);
3903 faces[0] =
Vec4I(0, 1, 2, 3);
3904 faces[1] =
Vec4I(7, 6, 5, 4);
3905 faces[2] =
Vec4I(4, 5, 1, 0);
3906 faces[3] =
Vec4I(6, 7, 3, 2);
3907 faces[4] =
Vec4I(0, 3, 7, 4);
3908 faces[5] =
Vec4I(1, 5, 6, 2);
3912 return meshToVolume<GridType>(mesh, xform, halfWidth, halfWidth);
3920 #endif // OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
void addLeaf(LeafNodeT *leaf)
Add the specified leaf to this tree, possibly creating a child branch in the process. If the leaf node already exists, replace it.
Definition: ValueAccessor.h:374
Index32 Index
Definition: Types.h:58
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:76
OPENVDB_API const Index32 INVALID_IDX
LeafNodeT * probeLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains voxel (x, y, z), or NULL if no such node exists...
Definition: ValueAccessor.h:424
NodeType * getNode()
Return the cached node of type NodeType. [Mainly for internal use].
Definition: ValueAccessor.h:349
Base class for tree-traversal iterators over tile and voxel values.
Definition: TreeIterator.h:658
bool normalize(T eps=T(1.0e-7))
this = normalized this
Definition: Vec3.h:348
void clear()
Remove all tiles from this tree and all nodes other than the root node.
Definition: Tree.h:610
#define OPENVDB_LOG_DEBUG(message)
In debug builds only, log a debugging message of the form 'someVar << "text" << ...'.
Definition: logging.h:45
LeafIter beginLeaf()
Return an iterator over all leaf nodes in this tree.
Definition: Tree.h:1143
Convert polygonal meshes that consist of quads and/or triangles into signed or unsigned distance fiel...
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:407
void clearAllAccessors()
Clear all registered accessors.
Definition: Tree.h:1477
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
Efficient multi-threaded replacement of the background values in tree.
Vec2< T > maxComponent(const Vec2< T > &v1, const Vec2< T > &v2)
Return component-wise maximum of the two vectors.
Definition: Vec2.h:520
uint32_t Index32
Definition: Types.h:56
OPENVDB_API Vec3d closestPointOnTriangleToPoint(const Vec3d &a, const Vec3d &b, const Vec3d &c, const Vec3d &p, Vec3d &uvw)
Closest Point on Triangle to Point. Given a triangle abc and a point p, return the point on abc close...
bool operator<(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition: Tuple.h:158
bool operator>(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition: Tuple.h:170
void setValueOn(const Coord &xyz, const ValueType &value)
Set the value of the voxel at the given coordinates and mark the voxel as active. ...
Definition: ValueAccessor.h:292
Vec3< T > cross(const Vec3< T > &v) const
Return the cross product of "this" vector and v;.
Definition: Vec3.h:232
Defined various multi-threaded utility functions for trees.
const Vec3T & min() const
Return a const reference to the minimum point of the BBox.
Definition: BBox.h:81
LeafNodeType * probeLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains voxel (x, y, z). If no such node exists, return NULL.
Definition: Tree.h:1664
Vec2< T > minComponent(const Vec2< T > &v1, const Vec2< T > &v2)
Return component-wise minimum of the two vectors.
Definition: Vec2.h:511
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: ValueAccessor.h:256
Tree< typename RootNodeType::template ValueConverter< Int32 >::Type > Type
Definition: Tree.h:223
const Vec3T & max() const
Return a const reference to the maximum point of the BBox.
Definition: BBox.h:84
Real GudonovsNormSqrd(bool isOutside, Real dP_xm, Real dP_xp, Real dP_ym, Real dP_yp, Real dP_zm, Real dP_zp)
Definition: FiniteDifference.h:353
#define OPENVDB_VERSION_NAME
Definition: version.h:43
bool isValueOn(const Coord &xyz) const
Return the active state of the voxel at the given coordinates.
Definition: ValueAccessor.h:263
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:52
void merge(Tree &other, MergePolicy=MERGE_ACTIVE_STATES)
Efficiently merge another tree into this tree using one of several schemes.
Definition: Tree.h:1774
Propagates the sign of distance values from the active voxels in the narrow band to the inactive valu...
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:324
const LeafNodeT * probeConstLeaf(const Coord &xyz) const
Return a pointer to the leaf node that contains voxel (x, y, z), or NULL if no such node exists...
Definition: ValueAccessor.h:429
Definition: Exceptions.h:39
RootNodeType::LeafNodeType LeafNodeType
Definition: Tree.h:211
Vec3< double > Vec3d
Definition: Vec3.h:643
Axis
Definition: Math.h:838
virtual void clear()
Remove all nodes from this cache, then reinsert the root node.
Definition: ValueAccessor.h:438
Type Pow2(Type x)
Return .
Definition: Math.h:514
Axis-aligned bounding box.
Definition: BBox.h:47
LeafNodeT * touchLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains voxel (x, y, z). If no such node exists, create one, but preserve the values and active states of all voxels.
Definition: ValueAccessor.h:393
const ValueT & getValue() const
Return the tile or voxel value to which this iterator is currently pointing.
Definition: TreeIterator.h:734
_RootNodeType RootNodeType
Definition: Tree.h:209
int32_t Int32
Definition: Types.h:60
math::Vec4< Index32 > Vec4I
Definition: Types.h:90
OPENVDB_API const Coord COORD_OFFSETS[26]
coordinate offset table for neighboring voxels
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:709
T dot(const Vec3< T > &v) const
Dot product.
Definition: Vec3.h:203
Definition: Exceptions.h:87
Vec3< float > Vec3s
Definition: Vec3.h:642
bool probeValue(const Coord &xyz, ValueType &value) const
Return the active state of the voxel as well as its value.
Definition: ValueAccessor.h:266
static const Real LEVEL_SET_HALF_WIDTH
Definition: Types.h:212
Base class for tree-traversal iterators over all leaf nodes (but not leaf voxels) ...
Definition: TreeIterator.h:1228