16 #ifndef OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
17 #define OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
30 #include <tbb/blocked_range.h>
31 #include <tbb/enumerable_thread_specific.h>
32 #include <tbb/parallel_for.h>
33 #include <tbb/parallel_reduce.h>
34 #include <tbb/partitioner.h>
35 #include <tbb/task_group.h>
36 #include <tbb/task_scheduler_init.h>
38 #include <boost/mpl/at.hpp>
39 #include <boost/mpl/int.hpp>
40 #include <boost/mpl/size.hpp>
48 #include <type_traits>
113 template <
typename Gr
idType,
typename MeshDataAdapter>
114 inline typename GridType::Ptr
118 float exteriorBandWidth = 3.0f,
119 float interiorBandWidth = 3.0f,
139 template <
typename Gr
idType,
typename MeshDataAdapter,
typename Interrupter>
140 inline typename GridType::Ptr
142 Interrupter& interrupter,
145 float exteriorBandWidth = 3.0f,
146 float interiorBandWidth = 3.0f,
164 template<
typename Po
intType,
typename PolygonType>
168 const std::vector<PolygonType>& polygons)
169 : mPointArray(points.empty() ? nullptr : &points[0])
170 , mPointArraySize(points.size())
171 , mPolygonArray(polygons.empty() ? nullptr : &polygons[0])
172 , mPolygonArraySize(polygons.size())
177 const PolygonType* polygonArray,
size_t polygonArraySize)
178 : mPointArray(pointArray)
179 , mPointArraySize(pointArraySize)
180 , mPolygonArray(polygonArray)
181 , mPolygonArraySize(polygonArraySize)
190 return (PolygonType::size == 3 || mPolygonArray[n][3] ==
util::INVALID_IDX) ? 3 : 4;
196 const PointType& p = mPointArray[mPolygonArray[n][int(v)]];
197 pos[0] = double(p[0]);
198 pos[1] = double(p[1]);
199 pos[2] = double(p[2]);
203 PointType
const *
const mPointArray;
204 size_t const mPointArraySize;
205 PolygonType
const *
const mPolygonArray;
206 size_t const mPolygonArraySize;
234 template<
typename Gr
idType>
235 inline typename GridType::Ptr
237 const openvdb::math::Transform& xform,
238 const std::vector<Vec3s>& points,
239 const std::vector<Vec3I>& triangles,
243 template<
typename Gr
idType,
typename Interrupter>
244 inline typename GridType::Ptr
246 Interrupter& interrupter,
247 const openvdb::math::Transform& xform,
248 const std::vector<Vec3s>& points,
249 const std::vector<Vec3I>& triangles,
268 template<
typename Gr
idType>
269 inline typename GridType::Ptr
271 const openvdb::math::Transform& xform,
272 const std::vector<Vec3s>& points,
273 const std::vector<Vec4I>& quads,
277 template<
typename Gr
idType,
typename Interrupter>
278 inline typename GridType::Ptr
280 Interrupter& interrupter,
281 const openvdb::math::Transform& xform,
282 const std::vector<Vec3s>& points,
283 const std::vector<Vec4I>& quads,
303 template<
typename Gr
idType>
304 inline typename GridType::Ptr
306 const openvdb::math::Transform& xform,
307 const std::vector<Vec3s>& points,
308 const std::vector<Vec3I>& triangles,
309 const std::vector<Vec4I>& quads,
313 template<
typename Gr
idType,
typename Interrupter>
314 inline typename GridType::Ptr
316 Interrupter& interrupter,
317 const openvdb::math::Transform& xform,
318 const std::vector<Vec3s>& points,
319 const std::vector<Vec3I>& triangles,
320 const std::vector<Vec4I>& quads,
342 template<
typename Gr
idType>
343 inline typename GridType::Ptr
345 const openvdb::math::Transform& xform,
346 const std::vector<Vec3s>& points,
347 const std::vector<Vec3I>& triangles,
348 const std::vector<Vec4I>& quads,
353 template<
typename Gr
idType,
typename Interrupter>
354 inline typename GridType::Ptr
356 Interrupter& interrupter,
357 const openvdb::math::Transform& xform,
358 const std::vector<Vec3s>& points,
359 const std::vector<Vec3I>& triangles,
360 const std::vector<Vec4I>& quads,
379 template<
typename Gr
idType>
380 inline typename GridType::Ptr
382 const openvdb::math::Transform& xform,
383 const std::vector<Vec3s>& points,
384 const std::vector<Vec3I>& triangles,
385 const std::vector<Vec4I>& quads,
389 template<
typename Gr
idType,
typename Interrupter>
390 inline typename GridType::Ptr
392 Interrupter& interrupter,
393 const openvdb::math::Transform& xform,
394 const std::vector<Vec3s>& points,
395 const std::vector<Vec3I>& triangles,
396 const std::vector<Vec4I>& quads,
409 template<
typename Gr
idType,
typename VecType>
410 inline typename GridType::Ptr
412 const openvdb::math::Transform& xform,
425 template <
typename FloatTreeT>
443 : mXDist(dist), mYDist(dist), mZDist(dist)
486 void convert(
const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList);
492 std::vector<Vec3d>& points, std::vector<Index32>& primitives);
511 namespace mesh_to_volume_internal {
513 template<
typename Po
intType>
518 : mPointsIn(pointsIn), mPointsOut(pointsOut), mXform(&xform)
522 void operator()(
const tbb::blocked_range<size_t>& range)
const {
526 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
528 const PointType& wsP = mPointsIn[n];
529 pos[0] = double(wsP[0]);
530 pos[1] = double(wsP[1]);
531 pos[2] = double(wsP[2]);
533 pos = mXform->worldToIndex(pos);
535 PointType& isP = mPointsOut[n];
536 isP[0] =
typename PointType::value_type(pos[0]);
537 isP[1] =
typename PointType::value_type(pos[1]);
538 isP[2] =
typename PointType::value_type(pos[2]);
548 template<
typename ValueType>
551 static ValueType
epsilon() {
return ValueType(1e-7); }
559 template<
typename TreeType>
571 : mDistTree(&lhsDistTree)
572 , mIdxTree(&lhsIdxTree)
573 , mRhsDistNodes(rhsDistNodes)
574 , mRhsIdxNodes(rhsIdxNodes)
578 void operator()(
const tbb::blocked_range<size_t>& range)
const {
583 using DistValueType =
typename LeafNodeType::ValueType;
584 using IndexValueType =
typename Int32LeafNodeType::ValueType;
586 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
588 const Coord& origin = mRhsDistNodes[n]->origin();
593 DistValueType* lhsDistData = lhsDistNode->buffer().data();
594 IndexValueType* lhsIdxData = lhsIdxNode->buffer().data();
596 const DistValueType* rhsDistData = mRhsDistNodes[n]->buffer().data();
597 const IndexValueType* rhsIdxData = mRhsIdxNodes[n]->buffer().data();
600 for (
Index32 offset = 0; offset < LeafNodeType::SIZE; ++offset) {
604 const DistValueType& lhsValue = lhsDistData[offset];
605 const DistValueType& rhsValue = rhsDistData[offset];
607 if (rhsValue < lhsValue) {
608 lhsDistNode->setValueOn(offset, rhsValue);
609 lhsIdxNode->setValueOn(offset, rhsIdxData[offset]);
611 lhsIdxNode->setValueOn(offset,
612 std::min(lhsIdxData[offset], rhsIdxData[offset]));
617 delete mRhsDistNodes[n];
618 delete mRhsIdxNodes[n];
624 TreeType *
const mDistTree;
625 Int32TreeType *
const mIdxTree;
627 LeafNodeType **
const mRhsDistNodes;
628 Int32LeafNodeType **
const mRhsIdxNodes;
635 template<
typename TreeType>
641 : mNodes(nodes.empty() ? nullptr : &nodes[0]), mCoordinates(coordinates)
645 void operator()(
const tbb::blocked_range<size_t>& range)
const {
646 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
647 Coord& origin = const_cast<Coord&>(mNodes[n]->origin());
648 mCoordinates[n] = origin;
649 origin[0] = static_cast<int>(n);
658 template<
typename TreeType>
664 : mNodes(nodes.empty() ? nullptr : &nodes[0]), mCoordinates(coordinates)
668 void operator()(
const tbb::blocked_range<size_t>& range)
const {
669 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
670 Coord& origin = const_cast<Coord&>(mNodes[n]->origin());
671 origin[0] = mCoordinates[n][0];
680 template<
typename TreeType>
687 size_t* offsets,
size_t numNodes,
const CoordBBox& bbox)
689 , mCoordinates(coordinates)
691 , mNumNodes(numNodes)
701 void operator()(
const tbb::blocked_range<size_t>& range)
const {
703 size_t* offsetsNextX = mOffsets;
704 size_t* offsetsPrevX = mOffsets + mNumNodes;
705 size_t* offsetsNextY = mOffsets + mNumNodes * 2;
706 size_t* offsetsPrevY = mOffsets + mNumNodes * 3;
707 size_t* offsetsNextZ = mOffsets + mNumNodes * 4;
708 size_t* offsetsPrevZ = mOffsets + mNumNodes * 5;
713 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
714 const Coord& origin = mCoordinates[n];
715 offsetsNextX[n] = findNeighbourNode(acc, origin,
Coord(LeafNodeType::DIM, 0, 0));
716 offsetsPrevX[n] = findNeighbourNode(acc, origin,
Coord(-LeafNodeType::DIM, 0, 0));
717 offsetsNextY[n] = findNeighbourNode(acc, origin,
Coord(0, LeafNodeType::DIM, 0));
718 offsetsPrevY[n] = findNeighbourNode(acc, origin,
Coord(0, -LeafNodeType::DIM, 0));
719 offsetsNextZ[n] = findNeighbourNode(acc, origin,
Coord(0, 0, LeafNodeType::DIM));
720 offsetsPrevZ[n] = findNeighbourNode(acc, origin,
Coord(0, 0, -LeafNodeType::DIM));
727 Coord ijk = start + step;
732 if (node)
return static_cast<size_t>(node->origin()[0]);
741 TreeType
const *
const mTree;
742 Coord const *
const mCoordinates;
743 size_t *
const mOffsets;
745 const size_t mNumNodes;
750 template<
typename TreeType>
759 mLeafNodes.reserve(tree.leafCount());
760 tree.getNodes(mLeafNodes);
762 if (mLeafNodes.empty())
return;
765 tree.evalLeafBoundingBox(bbox);
767 const tbb::blocked_range<size_t> range(0, mLeafNodes.size());
771 std::unique_ptr<Coord[]> coordinates{
new Coord[mLeafNodes.size()]};
772 tbb::parallel_for(range,
776 mOffsets.reset(
new size_t[mLeafNodes.size() * 6]);
780 tree, coordinates.get(), mOffsets.get(), mLeafNodes.size(), bbox));
786 size_t size()
const {
return mLeafNodes.size(); }
788 std::vector<LeafNodeType*>&
nodes() {
return mLeafNodes; }
789 const std::vector<LeafNodeType*>&
nodes()
const {
return mLeafNodes; }
793 const size_t*
offsetsPrevX()
const {
return mOffsets.get() + mLeafNodes.size(); }
795 const size_t*
offsetsNextY()
const {
return mOffsets.get() + mLeafNodes.size() * 2; }
796 const size_t*
offsetsPrevY()
const {
return mOffsets.get() + mLeafNodes.size() * 3; }
798 const size_t*
offsetsNextZ()
const {
return mOffsets.get() + mLeafNodes.size() * 4; }
799 const size_t*
offsetsPrevZ()
const {
return mOffsets.get() + mLeafNodes.size() * 5; }
802 std::vector<LeafNodeType*> mLeafNodes;
803 std::unique_ptr<size_t[]> mOffsets;
807 template<
typename TreeType>
820 : mStartNodeIndices(startNodeIndices.empty() ? nullptr : &startNodeIndices[0])
821 , mConnectivity(&connectivity)
826 void operator()(
const tbb::blocked_range<size_t>& range)
const {
828 std::vector<LeafNodeType*>& nodes = mConnectivity->nodes();
831 size_t idxA = 0, idxB = 1;
834 const size_t* nextOffsets = mConnectivity->offsetsNextZ();
835 const size_t* prevOffsets = mConnectivity->offsetsPrevZ();
841 step = LeafNodeType::DIM;
843 nextOffsets = mConnectivity->offsetsNextY();
844 prevOffsets = mConnectivity->offsetsPrevY();
846 }
else if (mAxis ==
X_AXIS) {
850 step = LeafNodeType::DIM * LeafNodeType::DIM;
852 nextOffsets = mConnectivity->offsetsNextX();
853 prevOffsets = mConnectivity->offsetsPrevX();
861 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
863 size_t startOffset = mStartNodeIndices[n];
864 size_t lastOffset = startOffset;
868 for (a = 0; a < int(LeafNodeType::DIM); ++a) {
869 for (b = 0; b < int(LeafNodeType::DIM); ++b) {
871 pos = LeafNodeType::coordToOffset(ijk);
872 size_t offset = startOffset;
875 while ( offset != ConnectivityTable::INVALID_OFFSET &&
876 traceVoxelLine(*nodes[offset], pos, step) ) {
879 offset = nextOffsets[offset];
884 while (offset != ConnectivityTable::INVALID_OFFSET) {
886 offset = nextOffsets[offset];
891 pos += step * (LeafNodeType::DIM - 1);
892 while ( offset != ConnectivityTable::INVALID_OFFSET &&
893 traceVoxelLine(*nodes[offset], pos, -step)) {
894 offset = prevOffsets[offset];
906 bool isOutside =
true;
908 for (
Index i = 0; i < LeafNodeType::DIM; ++i) {
916 if (!(dist >
ValueType(0.75))) isOutside =
false;
929 size_t const *
const mStartNodeIndices;
930 ConnectivityTable *
const mConnectivity;
936 template<
typename LeafNodeType>
940 using ValueType =
typename LeafNodeType::ValueType;
941 using Queue = std::deque<Index>;
944 ValueType* data = node.buffer().data();
948 for (
Index pos = 0; pos < LeafNodeType::SIZE; ++pos) {
949 if (data[pos] < 0.0) seedPoints.push_back(pos);
952 if (seedPoints.empty())
return;
955 for (Queue::iterator it = seedPoints.begin(); it != seedPoints.end(); ++it) {
956 ValueType& dist = data[*it];
963 Index pos(0), nextPos(0);
965 while (!seedPoints.empty()) {
967 pos = seedPoints.back();
968 seedPoints.pop_back();
970 ValueType& dist = data[pos];
972 if (!(dist < ValueType(0.0))) {
976 ijk = LeafNodeType::offsetToLocalCoord(pos);
979 nextPos = pos - LeafNodeType::DIM * LeafNodeType::DIM;
980 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
983 if (ijk[0] != (LeafNodeType::DIM - 1)) {
984 nextPos = pos + LeafNodeType::DIM * LeafNodeType::DIM;
985 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
989 nextPos = pos - LeafNodeType::DIM;
990 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
993 if (ijk[1] != (LeafNodeType::DIM - 1)) {
994 nextPos = pos + LeafNodeType::DIM;
995 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1000 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1003 if (ijk[2] != (LeafNodeType::DIM - 1)) {
1005 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1012 template<
typename LeafNodeType>
1016 bool updatedNode =
false;
1018 using ValueType =
typename LeafNodeType::ValueType;
1019 ValueType* data = node.buffer().data();
1023 bool updatedSign =
true;
1024 while (updatedSign) {
1026 updatedSign =
false;
1028 for (
Index pos = 0; pos < LeafNodeType::SIZE; ++pos) {
1030 ValueType& dist = data[pos];
1032 if (!(dist < ValueType(0.0)) && dist > ValueType(0.75)) {
1034 ijk = LeafNodeType::offsetToLocalCoord(pos);
1037 if (ijk[2] != 0 && data[pos - 1] < ValueType(0.0)) {
1039 dist = ValueType(-dist);
1042 }
else if (ijk[2] != (LeafNodeType::DIM - 1) && data[pos + 1] < ValueType(0.0)) {
1044 dist = ValueType(-dist);
1047 }
else if (ijk[1] != 0 && data[pos - LeafNodeType::DIM] < ValueType(0.0)) {
1049 dist = ValueType(-dist);
1052 }
else if (ijk[1] != (LeafNodeType::DIM - 1)
1053 && data[pos + LeafNodeType::DIM] < ValueType(0.0))
1056 dist = ValueType(-dist);
1059 }
else if (ijk[0] != 0
1060 && data[pos - LeafNodeType::DIM * LeafNodeType::DIM] < ValueType(0.0))
1063 dist = ValueType(-dist);
1066 }
else if (ijk[0] != (LeafNodeType::DIM - 1)
1067 && data[pos + LeafNodeType::DIM * LeafNodeType::DIM] < ValueType(0.0))
1070 dist = ValueType(-dist);
1075 updatedNode |= updatedSign;
1082 template<
typename TreeType>
1090 : mNodes(nodes.empty() ? nullptr : &nodes[0])
1091 , mChangedNodeMask(changedNodeMask)
1096 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1097 if (mChangedNodeMask[n]) {
1099 mChangedNodeMask[n] =
scanFill(*mNodes[n]);
1109 template<
typename ValueType>
1112 FillArray(ValueType* array,
const ValueType v) : mArray(array), mValue(v) { }
1115 const ValueType v = mValue;
1116 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1126 template<
typename ValueType>
1128 fillArray(ValueType* array,
const ValueType val,
const size_t length)
1130 const auto grainSize = std::max<size_t>(
1131 length / tbb::task_scheduler_init::default_num_threads(), 1024);
1132 const tbb::blocked_range<size_t> range(0, length, grainSize);
1137 template<
typename TreeType>
1145 const bool* changedNodeMask,
bool* changedVoxelMask)
1146 : mNodes(nodes.empty() ? nullptr : &nodes[0])
1147 , mChangedNodeMask(changedNodeMask)
1148 , mChangedVoxelMask(changedVoxelMask)
1153 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1155 if (mChangedNodeMask[n]) {
1156 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1158 ValueType* data = mNodes[n]->buffer().data();
1160 for (
Index pos = 0; pos < LeafNodeType::SIZE; ++pos) {
1176 template<
typename TreeType>
1185 bool* changedNodeMask,
bool* nodeMask,
bool* changedVoxelMask)
1186 : mConnectivity(&connectivity)
1187 , mChangedNodeMask(changedNodeMask)
1188 , mNodeMask(nodeMask)
1189 , mChangedVoxelMask(changedVoxelMask)
1195 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1197 if (!mChangedNodeMask[n]) {
1199 bool changedValue =
false;
1201 changedValue |= processZ(n,
true);
1202 changedValue |= processZ(n,
false);
1204 changedValue |= processY(n,
true);
1205 changedValue |= processY(n,
false);
1207 changedValue |= processX(n,
true);
1208 changedValue |= processX(n,
false);
1210 mNodeMask[n] = changedValue;
1218 const size_t offset =
1219 firstFace ? mConnectivity->offsetsPrevZ()[n] : mConnectivity->offsetsNextZ()[n];
1220 if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1222 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1224 const ValueType* lhsData = mConnectivity->nodes()[n]->buffer().data();
1225 const ValueType* rhsData = mConnectivity->nodes()[offset]->buffer().data();
1227 const Index lastOffset = LeafNodeType::DIM - 1;
1228 const Index lhsOffset =
1229 firstFace ? 0 : lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1231 Index tmpPos(0), pos(0);
1232 bool changedValue =
false;
1234 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
1235 tmpPos = x << (2 * LeafNodeType::LOG2DIM);
1236 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
1237 pos = tmpPos + (y << LeafNodeType::LOG2DIM);
1239 if (lhsData[pos + lhsOffset] >
ValueType(0.75)) {
1240 if (rhsData[pos + rhsOffset] <
ValueType(0.0)) {
1241 changedValue =
true;
1242 mask[pos + lhsOffset] =
true;
1248 return changedValue;
1256 const size_t offset =
1257 firstFace ? mConnectivity->offsetsPrevY()[n] : mConnectivity->offsetsNextY()[n];
1258 if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1260 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1262 const ValueType* lhsData = mConnectivity->nodes()[n]->buffer().data();
1263 const ValueType* rhsData = mConnectivity->nodes()[offset]->buffer().data();
1265 const Index lastOffset = LeafNodeType::DIM * (LeafNodeType::DIM - 1);
1266 const Index lhsOffset =
1267 firstFace ? 0 : lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1269 Index tmpPos(0), pos(0);
1270 bool changedValue =
false;
1272 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
1273 tmpPos = x << (2 * LeafNodeType::LOG2DIM);
1274 for (
Index z = 0; z < LeafNodeType::DIM; ++z) {
1277 if (lhsData[pos + lhsOffset] >
ValueType(0.75)) {
1278 if (rhsData[pos + rhsOffset] <
ValueType(0.0)) {
1279 changedValue =
true;
1280 mask[pos + lhsOffset] =
true;
1286 return changedValue;
1294 const size_t offset =
1295 firstFace ? mConnectivity->offsetsPrevX()[n] : mConnectivity->offsetsNextX()[n];
1296 if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1298 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1300 const ValueType* lhsData = mConnectivity->nodes()[n]->buffer().data();
1301 const ValueType* rhsData = mConnectivity->nodes()[offset]->buffer().data();
1303 const Index lastOffset = LeafNodeType::DIM * LeafNodeType::DIM * (LeafNodeType::DIM-1);
1304 const Index lhsOffset =
1305 firstFace ? 0 : lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1307 Index tmpPos(0), pos(0);
1308 bool changedValue =
false;
1310 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
1311 tmpPos = y << LeafNodeType::LOG2DIM;
1312 for (
Index z = 0; z < LeafNodeType::DIM; ++z) {
1315 if (lhsData[pos + lhsOffset] >
ValueType(0.75)) {
1316 if (rhsData[pos + rhsOffset] <
ValueType(0.0)) {
1317 changedValue =
true;
1318 mask[pos + lhsOffset] =
true;
1324 return changedValue;
1339 template<
typename TreeType,
typename MeshDataAdapter>
1353 std::vector<LeafNodeType*>& distNodes,
1354 const TreeType& distTree,
1357 : mDistNodes(distNodes.empty() ? nullptr : &distNodes[0])
1358 , mDistTree(&distTree)
1359 , mIndexTree(&indexTree)
1373 Index xPos(0), yPos(0);
1374 Coord ijk, nijk, nodeMin, nodeMax;
1375 Vec3d cp, xyz, nxyz, dir1, dir2;
1377 LocalData& localData = mLocalDataTable->local();
1380 if (!points) points.reset(
new Vec3d[LeafNodeType::SIZE * 2]);
1383 if (!mask) mask.reset(
new bool[LeafNodeType::SIZE]);
1386 typename LeafNodeType::ValueOnCIter it;
1388 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1394 const Int32* idxData = idxNode->buffer().data();
1396 nodeMin = node.origin();
1397 nodeMax = nodeMin.
offsetBy(LeafNodeType::DIM - 1);
1400 memset(mask.get(), 0,
sizeof(bool) * LeafNodeType::SIZE);
1402 for (it = node.cbeginValueOn(); it; ++it) {
1403 Index pos = it.pos();
1406 if (dist < 0.0 || dist > 0.75)
continue;
1408 ijk = node.offsetToGlobalCoord(pos);
1410 xyz[0] = double(ijk[0]);
1411 xyz[1] = double(ijk[1]);
1412 xyz[2] = double(ijk[2]);
1418 bool flipSign =
false;
1420 for (nijk[0] = bbox.
min()[0]; nijk[0] <= bbox.
max()[0] && !flipSign; ++nijk[0]) {
1421 xPos = (nijk[0] & (LeafNodeType::DIM - 1u)) << (2 * LeafNodeType::LOG2DIM);
1422 for (nijk[1]=bbox.
min()[1]; nijk[1] <= bbox.
max()[1] && !flipSign; ++nijk[1]) {
1423 yPos = xPos + ((nijk[1] & (LeafNodeType::DIM-1u)) << LeafNodeType::LOG2DIM);
1424 for (nijk[2] = bbox.
min()[2]; nijk[2] <= bbox.
max()[2]; ++nijk[2]) {
1425 pos = yPos + (nijk[2] & (LeafNodeType::DIM - 1u));
1427 const Int32& polyIdx = idxData[pos];
1432 const Index pointIndex = pos * 2;
1438 nxyz[0] = double(nijk[0]);
1439 nxyz[1] = double(nijk[1]);
1440 nxyz[2] = double(nijk[2]);
1442 Vec3d& point = points[pointIndex];
1444 point = closestPoint(nxyz, polyIdx);
1446 Vec3d& direction = points[pointIndex + 1];
1447 direction = nxyz - point;
1448 direction.normalize();
1451 dir1 = xyz - points[pointIndex];
1454 if (points[pointIndex + 1].dot(dir1) > 0.0) {
1465 for (
Int32 m = 0; m < 26; ++m) {
1469 nxyz[0] = double(nijk[0]);
1470 nxyz[1] = double(nijk[1]);
1471 nxyz[2] = double(nijk[2]);
1473 cp = closestPoint(nxyz, idxAcc.
getValue(nijk));
1481 if (dir2.dot(dir1) > 0.0) {
1497 Vec3d a, b, c, cp, uvw;
1499 const size_t polygon = size_t(polyIdx);
1500 mMesh->getIndexSpacePoint(polygon, 0, a);
1501 mMesh->getIndexSpacePoint(polygon, 1, b);
1502 mMesh->getIndexSpacePoint(polygon, 2, c);
1506 if (4 == mMesh->vertexCount(polygon)) {
1508 mMesh->getIndexSpacePoint(polygon, 3, b);
1512 if ((center - c).lengthSqr() < (center - cp).lengthSqr()) {
1521 LeafNodeType **
const mDistNodes;
1522 TreeType
const *
const mDistTree;
1523 Int32TreeType
const *
const mIndexTree;
1533 template<
typename LeafNodeType>
1537 using NodeT = LeafNodeType;
1539 const Coord ijk = NodeT::offsetToLocalCoord(pos);
1543 mask[0] = ijk[0] != (NodeT::DIM - 1);
1545 mask[1] = ijk[0] != 0;
1547 mask[2] = ijk[1] != (NodeT::DIM - 1);
1549 mask[3] = ijk[1] != 0;
1551 mask[4] = ijk[2] != (NodeT::DIM - 1);
1553 mask[5] = ijk[2] != 0;
1557 mask[6] = mask[0] && mask[5];
1559 mask[7] = mask[1] && mask[5];
1561 mask[8] = mask[0] && mask[4];
1563 mask[9] = mask[1] && mask[4];
1565 mask[10] = mask[0] && mask[2];
1567 mask[11] = mask[1] && mask[2];
1569 mask[12] = mask[0] && mask[3];
1571 mask[13] = mask[1] && mask[3];
1573 mask[14] = mask[3] && mask[4];
1575 mask[15] = mask[3] && mask[5];
1577 mask[16] = mask[2] && mask[4];
1579 mask[17] = mask[2] && mask[5];
1583 mask[18] = mask[1] && mask[3] && mask[5];
1585 mask[19] = mask[1] && mask[3] && mask[4];
1587 mask[20] = mask[0] && mask[3] && mask[4];
1589 mask[21] = mask[0] && mask[3] && mask[5];
1591 mask[22] = mask[1] && mask[2] && mask[5];
1593 mask[23] = mask[1] && mask[2] && mask[4];
1595 mask[24] = mask[0] && mask[2] && mask[4];
1597 mask[25] = mask[0] && mask[2] && mask[5];
1601 template<
typename Compare,
typename LeafNodeType>
1605 using NodeT = LeafNodeType;
1608 if (mask[5] && Compare::check(data[pos - 1]))
return true;
1610 if (mask[4] && Compare::check(data[pos + 1]))
return true;
1612 if (mask[3] && Compare::check(data[pos - NodeT::DIM]))
return true;
1614 if (mask[2] && Compare::check(data[pos + NodeT::DIM]))
return true;
1616 if (mask[1] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM]))
return true;
1618 if (mask[0] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM]))
return true;
1620 if (mask[6] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM]))
return true;
1622 if (mask[7] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - 1]))
return true;
1624 if (mask[8] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + 1]))
return true;
1626 if (mask[9] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + 1]))
return true;
1628 if (mask[10] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM]))
return true;
1630 if (mask[11] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM]))
return true;
1632 if (mask[12] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM]))
return true;
1634 if (mask[13] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM]))
return true;
1636 if (mask[14] && Compare::check(data[pos - NodeT::DIM + 1]))
return true;
1638 if (mask[15] && Compare::check(data[pos - NodeT::DIM - 1]))
return true;
1640 if (mask[16] && Compare::check(data[pos + NodeT::DIM + 1]))
return true;
1642 if (mask[17] && Compare::check(data[pos + NodeT::DIM - 1]))
return true;
1644 if (mask[18] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM - 1]))
return true;
1646 if (mask[19] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM + 1]))
return true;
1648 if (mask[20] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM + 1]))
return true;
1650 if (mask[21] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM - 1]))
return true;
1652 if (mask[22] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM - 1]))
return true;
1654 if (mask[23] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM + 1]))
return true;
1656 if (mask[24] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM + 1]))
return true;
1658 if (mask[25] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM - 1]))
return true;
1664 template<
typename Compare,
typename AccessorType>
1668 for (
Int32 m = 0; m < 26; ++m) {
1678 template<
typename TreeType>
1688 , mNodes(nodes.empty() ? nullptr : &nodes[0])
1695 bool neighbourMask[26];
1697 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1702 typename LeafNodeType::ValueOnCIter it;
1703 for (it = node.cbeginValueOn(); it; ++it) {
1705 const Index pos = it.pos();
1708 if (dist < 0.0 || dist > 0.75)
continue;
1711 maskNodeInternalNeighbours<LeafNodeType>(pos, neighbourMask);
1713 const bool hasNegativeNeighbour =
1714 checkNeighbours<IsNegative, LeafNodeType>(pos, data, neighbourMask) ||
1715 checkNeighbours<IsNegative>(node.offsetToGlobalCoord(pos), acc, neighbourMask);
1717 if (!hasNegativeNeighbour) {
1730 template<
typename TreeType>
1741 : mNodes(nodes.empty() ? nullptr : &nodes[0])
1742 , mDistTree(&distTree)
1743 , mIndexTree(&indexTree)
1751 bool neighbourMask[26];
1753 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1756 ValueType* data = distNode.buffer().data();
1758 typename Int32TreeType::LeafNodeType* idxNode =
1761 typename LeafNodeType::ValueOnCIter it;
1762 for (it = distNode.cbeginValueOn(); it; ++it) {
1764 const Index pos = it.pos();
1766 if (!(data[pos] > 0.75))
continue;
1769 maskNodeInternalNeighbours<LeafNodeType>(pos, neighbourMask);
1771 const bool hasBoundaryNeighbour =
1772 checkNeighbours<Comp, LeafNodeType>(pos, data, neighbourMask) ||
1773 checkNeighbours<Comp>(distNode.offsetToGlobalCoord(pos),distAcc,neighbourMask);
1775 if (!hasBoundaryNeighbour) {
1776 distNode.setValueOff(pos);
1777 idxNode->setValueOff(pos);
1792 template<
typename NodeType>
1799 using NodeMaskType =
typename NodeType::NodeMaskType;
1801 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1802 const_cast<NodeMaskType&>(mNodes[n]->getChildMask()).setOff();
1810 template<
typename TreeType>
1814 using RootNodeType =
typename TreeType::RootNodeType;
1815 using NodeChainType =
typename RootNodeType::NodeChainType;
1816 using InternalNodeType =
typename boost::mpl::at<NodeChainType, boost::mpl::int_<1> >::type;
1818 std::vector<InternalNodeType*> nodes;
1819 tree.getNodes(nodes);
1821 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
1826 template<
typename TreeType>
1832 std::vector<LeafNodeType*>& overlappingNodes)
1833 : mLhsTree(&lhsTree)
1834 , mRhsTree(&rhsTree)
1835 , mNodes(&overlappingNodes)
1841 std::vector<LeafNodeType*> rhsLeafNodes;
1843 rhsLeafNodes.reserve(mRhsTree->leafCount());
1846 mRhsTree->stealNodes(rhsLeafNodes);
1850 for (
size_t n = 0, N = rhsLeafNodes.size(); n < N; ++n) {
1851 if (!acc.
probeLeaf(rhsLeafNodes[n]->origin())) {
1854 mNodes->push_back(rhsLeafNodes[n]);
1860 TreeType *
const mLhsTree;
1861 TreeType *
const mRhsTree;
1862 std::vector<LeafNodeType*> *
const mNodes;
1866 template<
typename DistTreeType,
typename IndexTreeType>
1869 DistTreeType& rhsDist, IndexTreeType& rhsIdx)
1871 using DistLeafNodeType =
typename DistTreeType::LeafNodeType;
1872 using IndexLeafNodeType =
typename IndexTreeType::LeafNodeType;
1874 std::vector<DistLeafNodeType*> overlappingDistNodes;
1875 std::vector<IndexLeafNodeType*> overlappingIdxNodes;
1878 tbb::task_group tasks;
1884 if (!overlappingDistNodes.empty() && !overlappingIdxNodes.empty()) {
1885 tbb::parallel_for(tbb::blocked_range<size_t>(0, overlappingDistNodes.size()),
1887 &overlappingDistNodes[0], &overlappingIdxNodes[0]));
1897 template<
typename TreeType>
1900 using Ptr = std::unique_ptr<VoxelizationData>;
1904 using UCharTreeType =
typename TreeType::template ValueConverter<unsigned char>::Type;
1915 , indexAcc(indexTree)
1916 , primIdTree(MaxPrimId)
1917 , primIdAcc(primIdTree)
1947 if (mPrimCount == MaxPrimId || primIdTree.leafCount() > 1000) {
1949 primIdTree.root().clear();
1950 primIdTree.clearAllAccessors();
1951 assert(mPrimCount == 0);
1954 return mPrimCount++;
1959 enum { MaxPrimId = 100 };
1961 unsigned char mPrimCount;
1965 template<
typename TreeType,
typename MeshDataAdapter,
typename Interrupter = util::NullInterrupter>
1971 using DataTable = tbb::enumerable_thread_specific<typename VoxelizationDataType::Ptr>;
1975 Interrupter* interrupter =
nullptr)
1976 : mDataTable(&dataTable)
1978 , mInterrupter(interrupter)
1989 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1992 tbb::task::self().cancel_group_execution();
1996 const size_t numVerts = mMesh->vertexCount(n);
1999 if (numVerts == 3 || numVerts == 4) {
2001 prim.index =
Int32(n);
2003 mMesh->getIndexSpacePoint(n, 0, prim.a);
2004 mMesh->getIndexSpacePoint(n, 1, prim.b);
2005 mMesh->getIndexSpacePoint(n, 2, prim.c);
2007 evalTriangle(prim, *dataPtr);
2009 if (numVerts == 4) {
2010 mMesh->getIndexSpacePoint(n, 3, prim.b);
2011 evalTriangle(prim, *dataPtr);
2019 bool wasInterrupted()
const {
return mInterrupter && mInterrupter->wasInterrupted(); }
2021 struct Triangle {
Vec3d a, b, c;
Int32 index; };
2025 enum { POLYGON_LIMIT = 1000 };
2027 SubTask(
const Triangle& prim, DataTable& dataTable,
2028 int subdivisionCount,
size_t polygonCount,
2029 Interrupter* interrupter =
nullptr)
2030 : mLocalDataTable(&dataTable)
2032 , mSubdivisionCount(subdivisionCount)
2033 , mPolygonCount(polygonCount)
2034 , mInterrupter(interrupter)
2038 void operator()()
const
2040 if (mSubdivisionCount <= 0 || mPolygonCount >= POLYGON_LIMIT) {
2042 typename VoxelizationDataType::Ptr& dataPtr = mLocalDataTable->local();
2043 if (!dataPtr) dataPtr.reset(
new VoxelizationDataType());
2045 voxelizeTriangle(mPrim, *dataPtr, mInterrupter);
2047 }
else if (!(mInterrupter && mInterrupter->wasInterrupted())) {
2048 spawnTasks(mPrim, *mLocalDataTable, mSubdivisionCount, mPolygonCount, mInterrupter);
2052 DataTable *
const mLocalDataTable;
2053 Triangle
const mPrim;
2054 int const mSubdivisionCount;
2055 size_t const mPolygonCount;
2056 Interrupter *
const mInterrupter;
2059 inline static int evalSubdivisionCount(
const Triangle& prim)
2061 const double ax = prim.a[0], bx = prim.b[0], cx = prim.c[0];
2064 const double ay = prim.a[1], by = prim.b[1], cy = prim.c[1];
2067 const double az = prim.a[2], bz = prim.b[2], cz = prim.c[2];
2070 return int(
std::max(dx,
std::max(dy, dz)) /
double(TreeType::LeafNodeType::DIM * 2));
2073 void evalTriangle(
const Triangle& prim, VoxelizationDataType& data)
const
2075 const size_t polygonCount = mMesh->polygonCount();
2076 const int subdivisionCount =
2077 polygonCount < SubTask::POLYGON_LIMIT ? evalSubdivisionCount(prim) : 0;
2079 if (subdivisionCount <= 0) {
2080 voxelizeTriangle(prim, data, mInterrupter);
2082 spawnTasks(prim, *mDataTable, subdivisionCount, polygonCount, mInterrupter);
2086 static void spawnTasks(
2087 const Triangle& mainPrim,
2088 DataTable& dataTable,
2089 int subdivisionCount,
2090 size_t polygonCount,
2091 Interrupter*
const interrupter)
2093 subdivisionCount -= 1;
2096 tbb::task_group tasks;
2098 const Vec3d ac = (mainPrim.a + mainPrim.c) * 0.5;
2099 const Vec3d bc = (mainPrim.b + mainPrim.c) * 0.5;
2100 const Vec3d ab = (mainPrim.a + mainPrim.b) * 0.5;
2103 prim.index = mainPrim.index;
2105 prim.a = mainPrim.a;
2108 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount, interrupter));
2113 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount, interrupter));
2116 prim.b = mainPrim.b;
2118 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount, interrupter));
2122 prim.c = mainPrim.c;
2123 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount, interrupter));
2128 static void voxelizeTriangle(
const Triangle& prim, VoxelizationDataType& data, Interrupter*
const interrupter)
2130 std::deque<Coord> coordList;
2134 coordList.push_back(ijk);
2139 updateDistance(ijk, prim, data);
2141 unsigned char primId = data.getNewPrimId();
2142 data.primIdAcc.setValueOnly(ijk, primId);
2144 while (!coordList.empty()) {
2145 if (interrupter && interrupter->wasInterrupted()) {
2146 tbb::task::self().cancel_group_execution();
2149 for (
Int32 pass = 0; pass < 1048576 && !coordList.empty(); ++pass) {
2150 ijk = coordList.back();
2151 coordList.pop_back();
2153 for (
Int32 i = 0; i < 26; ++i) {
2155 if (primId != data.primIdAcc.getValue(nijk)) {
2156 data.primIdAcc.setValueOnly(nijk, primId);
2157 if(updateDistance(nijk, prim, data)) coordList.push_back(nijk);
2164 static bool updateDistance(
const Coord& ijk,
const Triangle& prim, VoxelizationDataType& data)
2166 Vec3d uvw, voxelCenter(ijk[0], ijk[1], ijk[2]);
2168 using ValueType =
typename TreeType::ValueType;
2170 const ValueType dist = ValueType((voxelCenter -
2175 if (std::isnan(dist))
2178 const ValueType oldDist = data.distAcc.getValue(ijk);
2180 if (dist < oldDist) {
2181 data.distAcc.setValue(ijk, dist);
2182 data.indexAcc.setValue(ijk, prim.index);
2186 data.indexAcc.setValueOnly(ijk,
std::min(prim.index, data.indexAcc.getValue(ijk)));
2189 return !(dist > 0.75);
2192 DataTable *
const mDataTable;
2194 Interrupter *
const mInterrupter;
2201 template<
typename TreeType>
2207 using BoolTreeType =
typename TreeType::template ValueConverter<bool>::Type;
2211 std::vector<BoolLeafNodeType*>& lhsNodes)
2212 : mRhsTree(&rhsTree), mLhsNodes(lhsNodes.empty() ? nullptr : &lhsNodes[0])
2220 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2225 if (rhsNode) lhsNode->topologyDifference(*rhsNode,
false);
2230 TreeType
const *
const mRhsTree;
2231 BoolLeafNodeType **
const mLhsNodes;
2235 template<
typename LeafNodeTypeA,
typename LeafNodeTypeB>
2239 : mNodesA(nodesA.empty() ? nullptr : &nodesA[0])
2240 , mNodesB(nodesB.empty() ? nullptr : &nodesB[0])
2245 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2246 mNodesA[n]->topologyUnion(*mNodesB[n]);
2251 LeafNodeTypeA **
const mNodesA;
2252 LeafNodeTypeB **
const mNodesB;
2256 template<
typename TreeType>
2261 using BoolTreeType =
typename TreeType::template ValueConverter<bool>::Type;
2265 std::vector<LeafNodeType*>& nodes)
2267 , mNodes(nodes.empty() ? nullptr : &nodes[0])
2268 , mLocalMaskTree(false)
2269 , mMaskTree(&maskTree)
2275 , mNodes(rhs.mNodes)
2276 , mLocalMaskTree(false)
2277 , mMaskTree(&mLocalMaskTree)
2283 using Iterator =
typename LeafNodeType::ValueOnCIter;
2288 Coord ijk, nijk, localCorod;
2291 for (
size_t n = range.begin(); n != range.end(); ++n) {
2295 CoordBBox bbox = node.getNodeBoundingBox();
2300 for (Iterator it = node.cbeginValueOn(); it; ++it) {
2301 ijk = it.getCoord();
2304 localCorod = LeafNodeType::offsetToLocalCoord(pos);
2306 if (localCorod[2] <
int(LeafNodeType::DIM - 1)) {
2308 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2314 if (localCorod[2] > 0) {
2316 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2322 if (localCorod[1] <
int(LeafNodeType::DIM - 1)) {
2323 npos = pos + LeafNodeType::DIM;
2324 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2330 if (localCorod[1] > 0) {
2331 npos = pos - LeafNodeType::DIM;
2332 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2338 if (localCorod[0] <
int(LeafNodeType::DIM - 1)) {
2339 npos = pos + LeafNodeType::DIM * LeafNodeType::DIM;
2340 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2346 if (localCorod[0] > 0) {
2347 npos = pos - LeafNodeType::DIM * LeafNodeType::DIM;
2348 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2360 TreeType
const *
const mTree;
2361 LeafNodeType **
const mNodes;
2363 BoolTreeType mLocalMaskTree;
2364 BoolTreeType *
const mMaskTree;
2369 template<
typename TreeType,
typename MeshDataAdapter>
2377 using BoolTreeType =
typename TreeType::template ValueConverter<bool>::Type;
2388 : idx(idx_), x(x_), y(y_), z(z_), dist(dist_)
2398 std::vector<BoolLeafNodeType*>& maskNodes,
2406 : mMaskNodes(maskNodes.empty() ? nullptr : &maskNodes[0])
2407 , mMaskTree(&maskTree)
2408 , mDistTree(&distTree)
2409 , mIndexTree(&indexTree)
2411 , mNewMaskTree(false)
2413 , mUpdatedDistNodes()
2415 , mUpdatedIndexNodes()
2416 , mExteriorBandWidth(exteriorBandWidth)
2417 , mInteriorBandWidth(interiorBandWidth)
2418 , mVoxelSize(voxelSize)
2423 : mMaskNodes(rhs.mMaskNodes)
2424 , mMaskTree(rhs.mMaskTree)
2425 , mDistTree(rhs.mDistTree)
2426 , mIndexTree(rhs.mIndexTree)
2428 , mNewMaskTree(false)
2430 , mUpdatedDistNodes()
2432 , mUpdatedIndexNodes()
2433 , mExteriorBandWidth(rhs.mExteriorBandWidth)
2434 , mInteriorBandWidth(rhs.mInteriorBandWidth)
2435 , mVoxelSize(rhs.mVoxelSize)
2441 mDistNodes.insert(mDistNodes.end(), rhs.mDistNodes.begin(), rhs.mDistNodes.end());
2442 mIndexNodes.insert(mIndexNodes.end(), rhs.mIndexNodes.begin(), rhs.mIndexNodes.end());
2444 mUpdatedDistNodes.insert(mUpdatedDistNodes.end(),
2445 rhs.mUpdatedDistNodes.begin(), rhs.mUpdatedDistNodes.end());
2447 mUpdatedIndexNodes.insert(mUpdatedIndexNodes.end(),
2448 rhs.mUpdatedIndexNodes.begin(), rhs.mUpdatedIndexNodes.end());
2450 mNewMaskTree.merge(rhs.mNewMaskTree);
2459 std::vector<Fragment> fragments;
2460 fragments.reserve(256);
2462 std::unique_ptr<LeafNodeType> newDistNodePt;
2463 std::unique_ptr<Int32LeafNodeType> newIndexNodePt;
2465 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2468 if (maskNode.isEmpty())
continue;
2472 const Coord& origin = maskNode.origin();
2477 assert(!distNodePt == !indexNodePt);
2479 bool usingNewNodes =
false;
2481 if (!distNodePt && !indexNodePt) {
2485 if (!newDistNodePt.get() && !newIndexNodePt.get()) {
2486 newDistNodePt.reset(
new LeafNodeType(origin, backgroundDist));
2490 if ((backgroundDist <
ValueType(0.0)) !=
2491 (newDistNodePt->getValue(0) <
ValueType(0.0))) {
2492 newDistNodePt->buffer().fill(backgroundDist);
2495 newDistNodePt->setOrigin(origin);
2496 newIndexNodePt->setOrigin(origin);
2499 distNodePt = newDistNodePt.get();
2500 indexNodePt = newIndexNodePt.get();
2502 usingNewNodes =
true;
2509 for (
typename BoolLeafNodeType::ValueOnIter it = maskNode.beginValueOn(); it; ++it) {
2510 bbox.
expand(it.getCoord());
2515 gatherFragments(fragments, bbox, distAcc, indexAcc);
2520 bbox = maskNode.getNodeBoundingBox();
2522 bool updatedLeafNodes =
false;
2524 for (
typename BoolLeafNodeType::ValueOnIter it = maskNode.beginValueOn(); it; ++it) {
2526 const Coord ijk = it.getCoord();
2528 if (updateVoxel(ijk, 5, fragments, *distNodePt, *indexNodePt, &updatedLeafNodes)) {
2530 for (
Int32 i = 0; i < 6; ++i) {
2533 mask.setOn(BoolLeafNodeType::coordToOffset(nijk));
2539 for (
Int32 i = 6; i < 26; ++i) {
2542 mask.setOn(BoolLeafNodeType::coordToOffset(nijk));
2548 if (updatedLeafNodes) {
2551 mask -= indexNodePt->getValueMask();
2553 for (
typename NodeMaskType::OnIterator it = mask.beginOn(); it; ++it) {
2555 const Index pos = it.pos();
2556 const Coord ijk = maskNode.origin() + LeafNodeType::offsetToLocalCoord(pos);
2558 if (updateVoxel(ijk, 6, fragments, *distNodePt, *indexNodePt)) {
2559 for (
Int32 i = 0; i < 6; ++i) {
2566 if (usingNewNodes) {
2567 newDistNodePt->topologyUnion(*newIndexNodePt);
2568 mDistNodes.push_back(newDistNodePt.release());
2569 mIndexNodes.push_back(newIndexNodePt.release());
2571 mUpdatedDistNodes.push_back(distNodePt);
2572 mUpdatedIndexNodes.push_back(indexNodePt);
2592 gatherFragments(std::vector<Fragment>& fragments,
const CoordBBox& bbox,
2596 const Coord nodeMin = bbox.
min() & ~(LeafNodeType::DIM - 1);
2597 const Coord nodeMax = bbox.
max() & ~(LeafNodeType::DIM - 1);
2602 for (ijk[0] = nodeMin[0]; ijk[0] <= nodeMax[0]; ijk[0] += LeafNodeType::DIM) {
2603 for (ijk[1] = nodeMin[1]; ijk[1] <= nodeMax[1]; ijk[1] += LeafNodeType::DIM) {
2604 for (ijk[2] = nodeMin[2]; ijk[2] <= nodeMax[2]; ijk[2] += LeafNodeType::DIM) {
2605 if (LeafNodeType* distleaf = distAcc.
probeLeaf(ijk)) {
2608 ijk.
offsetBy(LeafNodeType::DIM - 1));
2609 gatherFragments(fragments, region, *distleaf, *indexAcc.
probeLeaf(ijk));
2615 std::sort(fragments.begin(), fragments.end());
2619 gatherFragments(std::vector<Fragment>& fragments,
const CoordBBox& bbox,
2620 const LeafNodeType& distLeaf,
const Int32LeafNodeType& idxLeaf)
const
2622 const typename LeafNodeType::NodeMaskType& mask = distLeaf.getValueMask();
2623 const ValueType* distData = distLeaf.buffer().data();
2624 const Int32* idxData = idxLeaf.buffer().data();
2626 for (
int x = bbox.min()[0]; x <= bbox.max()[0]; ++x) {
2627 const Index xPos = (x & (LeafNodeType::DIM - 1u)) << (2 * LeafNodeType::LOG2DIM);
2628 for (
int y = bbox.min()[1]; y <= bbox.max()[1]; ++y) {
2629 const Index yPos = xPos + ((y & (LeafNodeType::DIM - 1u)) << LeafNodeType::LOG2DIM);
2630 for (
int z = bbox.min()[2]; z <= bbox.max()[2]; ++z) {
2631 const Index pos = yPos + (z & (LeafNodeType::DIM - 1u));
2632 if (mask.isOn(pos)) {
2633 fragments.push_back(Fragment(idxData[pos],x,y,z, std::abs(distData[pos])));
2643 computeDistance(
const Coord& ijk,
const Int32 manhattanLimit,
2644 const std::vector<Fragment>& fragments,
Int32& closestPrimIdx)
const
2646 Vec3d a, b, c, uvw, voxelCenter(ijk[0], ijk[1], ijk[2]);
2650 for (
size_t n = 0, N = fragments.size(); n < N; ++n) {
2652 const Fragment& fragment = fragments[n];
2653 if (lastIdx == fragment.idx)
continue;
2655 const Int32 dx = std::abs(fragment.x - ijk[0]);
2656 const Int32 dy = std::abs(fragment.y - ijk[1]);
2657 const Int32 dz = std::abs(fragment.z - ijk[2]);
2659 const Int32 manhattan = dx + dy + dz;
2660 if (manhattan > manhattanLimit)
continue;
2662 lastIdx = fragment.idx;
2664 const size_t polygon = size_t(lastIdx);
2666 mMesh->getIndexSpacePoint(polygon, 0, a);
2667 mMesh->getIndexSpacePoint(polygon, 1, b);
2668 mMesh->getIndexSpacePoint(polygon, 2, c);
2670 primDist = (voxelCenter -
2674 if (4 == mMesh->vertexCount(polygon)) {
2676 mMesh->getIndexSpacePoint(polygon, 3, b);
2679 a, b, c, voxelCenter, uvw)).lengthSqr();
2681 if (tmpDist < primDist) primDist = tmpDist;
2684 if (primDist < dist) {
2686 closestPrimIdx = lastIdx;
2690 return ValueType(std::sqrt(dist)) * mVoxelSize;
2696 updateVoxel(
const Coord& ijk,
const Int32 manhattanLimit,
2697 const std::vector<Fragment>& fragments,
2698 LeafNodeType& distLeaf, Int32LeafNodeType& idxLeaf,
bool* updatedLeafNodes =
nullptr)
2700 Int32 closestPrimIdx = 0;
2701 const ValueType distance = computeDistance(ijk, manhattanLimit, fragments, closestPrimIdx);
2703 const Index pos = LeafNodeType::coordToOffset(ijk);
2704 const bool inside = distLeaf.getValue(pos) < ValueType(0.0);
2706 bool activateNeighbourVoxels =
false;
2708 if (!inside && distance < mExteriorBandWidth) {
2709 if (updatedLeafNodes) *updatedLeafNodes =
true;
2710 activateNeighbourVoxels = (distance + mVoxelSize) < mExteriorBandWidth;
2711 distLeaf.setValueOnly(pos, distance);
2712 idxLeaf.setValueOn(pos, closestPrimIdx);
2713 }
else if (inside && distance < mInteriorBandWidth) {
2714 if (updatedLeafNodes) *updatedLeafNodes =
true;
2715 activateNeighbourVoxels = (distance + mVoxelSize) < mInteriorBandWidth;
2716 distLeaf.setValueOnly(pos, -distance);
2717 idxLeaf.setValueOn(pos, closestPrimIdx);
2720 return activateNeighbourVoxels;
2725 BoolLeafNodeType **
const mMaskNodes;
2726 BoolTreeType *
const mMaskTree;
2727 TreeType *
const mDistTree;
2728 Int32TreeType *
const mIndexTree;
2732 BoolTreeType mNewMaskTree;
2734 std::vector<LeafNodeType*> mDistNodes, mUpdatedDistNodes;
2735 std::vector<Int32LeafNodeType*> mIndexNodes, mUpdatedIndexNodes;
2737 const ValueType mExteriorBandWidth, mInteriorBandWidth, mVoxelSize;
2741 template<
typename TreeType>
2745 AddNodes(TreeType& tree, std::vector<LeafNodeType*>& nodes)
2746 : mTree(&tree) , mNodes(&nodes)
2752 std::vector<LeafNodeType*>& nodes = *mNodes;
2753 for (
size_t n = 0, N = nodes.size(); n < N; ++n) {
2763 template<
typename TreeType,
typename Int32TreeType,
typename BoolTreeType,
typename MeshDataAdapter>
2767 Int32TreeType& indexTree,
2768 BoolTreeType& maskTree,
2769 std::vector<typename BoolTreeType::LeafNodeType*>& maskNodes,
2771 typename TreeType::ValueType exteriorBandWidth,
2772 typename TreeType::ValueType interiorBandWidth,
2773 typename TreeType::ValueType voxelSize)
2776 distTree, indexTree, mesh, exteriorBandWidth, interiorBandWidth, voxelSize);
2778 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, maskNodes.size()), expandOp);
2780 tbb::parallel_for(tbb::blocked_range<size_t>(0, expandOp.
updatedIndexNodes().size()),
2784 tbb::task_group tasks;
2798 template<
typename TreeType>
2807 , mVoxelSize(voxelSize)
2808 , mUnsigned(unsignedDist)
2814 typename LeafNodeType::ValueOnIter iter;
2816 const bool udf = mUnsigned;
2817 const ValueType w[2] = { -mVoxelSize, mVoxelSize };
2819 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2821 for (iter = mNodes[n]->beginValueOn(); iter; ++iter) {
2822 ValueType& val = const_cast<ValueType&>(iter.getValue());
2823 val = w[udf || (val <
ValueType(0.0))] * std::sqrt(std::abs(val));
2829 LeafNodeType * *
const mNodes;
2830 const ValueType mVoxelSize;
2831 const bool mUnsigned;
2836 template<
typename TreeType>
2844 : mNodes(nodes.empty() ? nullptr : &nodes[0])
2845 , mExBandWidth(exBandWidth)
2846 , mInBandWidth(inBandWidth)
2852 typename LeafNodeType::ValueOnIter iter;
2856 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2858 for (iter = mNodes[n]->beginValueOn(); iter; ++iter) {
2860 ValueType& val = const_cast<ValueType&>(iter.getValue());
2862 const bool inside = val <
ValueType(0.0);
2864 if (inside && !(val > inVal)) {
2867 }
else if (!inside && !(val < exVal)) {
2876 LeafNodeType * *
const mNodes;
2877 const ValueType mExBandWidth, mInBandWidth;
2881 template<
typename TreeType>
2888 : mNodes(nodes.empty() ? nullptr : &nodes[0]), mOffset(offset)
2896 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2898 typename LeafNodeType::ValueOnIter iter = mNodes[n]->beginValueOn();
2900 for (; iter; ++iter) {
2901 ValueType& val = const_cast<ValueType&>(iter.getValue());
2908 LeafNodeType * *
const mNodes;
2909 const ValueType mOffset;
2913 template<
typename TreeType>
2919 Renormalize(
const TreeType& tree,
const std::vector<LeafNodeType*>& nodes,
2922 , mNodes(nodes.empty() ? nullptr : &nodes[0])
2924 , mVoxelSize(voxelSize)
2939 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2941 ValueType* bufferData = &mBuffer[n * LeafNodeType::SIZE];
2943 typename LeafNodeType::ValueOnCIter iter = mNodes[n]->cbeginValueOn();
2944 for (; iter; ++iter) {
2948 ijk = iter.getCoord();
2963 bufferData[iter.pos()] = phi0 - dx * S * diff;
2969 TreeType
const *
const mTree;
2970 LeafNodeType
const *
const *
const mNodes;
2971 ValueType *
const mBuffer;
2973 const ValueType mVoxelSize;
2977 template<
typename TreeType>
2984 : mNodes(nodes.empty() ? nullptr : &nodes[0]), mBuffer(buffer)
2990 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2992 const ValueType* bufferData = &mBuffer[n * LeafNodeType::SIZE];
2994 typename LeafNodeType::ValueOnIter iter = mNodes[n]->beginValueOn();
2996 for (; iter; ++iter) {
2997 ValueType& val = const_cast<ValueType&>(iter.getValue());
2998 val =
std::min(val, bufferData[iter.pos()]);
3004 LeafNodeType * *
const mNodes;
3005 ValueType
const *
const mBuffer;
3017 template <
typename FloatTreeT>
3023 ConnectivityTable nodeConnectivity(tree);
3025 std::vector<size_t> zStartNodes, yStartNodes, xStartNodes;
3027 for (
size_t n = 0; n < nodeConnectivity.size(); ++n) {
3028 if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevX()[n]) {
3029 xStartNodes.push_back(n);
3032 if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevY()[n]) {
3033 yStartNodes.push_back(n);
3036 if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevZ()[n]) {
3037 zStartNodes.push_back(n);
3043 tbb::parallel_for(tbb::blocked_range<size_t>(0, zStartNodes.size()),
3046 tbb::parallel_for(tbb::blocked_range<size_t>(0, yStartNodes.size()),
3049 tbb::parallel_for(tbb::blocked_range<size_t>(0, xStartNodes.size()),
3052 const size_t numLeafNodes = nodeConnectivity.size();
3053 const size_t numVoxels = numLeafNodes * FloatTreeT::LeafNodeType::SIZE;
3055 std::unique_ptr<bool[]> changedNodeMaskA{
new bool[numLeafNodes]};
3056 std::unique_ptr<bool[]> changedNodeMaskB{
new bool[numLeafNodes]};
3057 std::unique_ptr<bool[]> changedVoxelMask{
new bool[numVoxels]};
3063 const tbb::blocked_range<size_t> nodeRange(0, numLeafNodes);
3065 bool nodesUpdated =
false;
3068 nodeConnectivity.nodes(), changedNodeMaskA.get()));
3071 nodeConnectivity, changedNodeMaskA.get(), changedNodeMaskB.get(),
3072 changedVoxelMask.get()));
3074 changedNodeMaskA.swap(changedNodeMaskB);
3076 nodesUpdated =
false;
3077 for (
size_t n = 0; n < numLeafNodes; ++n) {
3078 nodesUpdated |= changedNodeMaskA[n];
3079 if (nodesUpdated)
break;
3084 nodeConnectivity.nodes(), changedNodeMaskA.get(), changedVoxelMask.get()));
3086 }
while (nodesUpdated);
3094 template <
typename Gr
idType,
typename MeshDataAdapter,
typename Interrupter>
3095 inline typename GridType::Ptr
3097 Interrupter& interrupter,
3100 float exteriorBandWidth,
3101 float interiorBandWidth,
3105 using GridTypePtr =
typename GridType::Ptr;
3106 using TreeType =
typename GridType::TreeType;
3107 using LeafNodeType =
typename TreeType::LeafNodeType;
3108 using ValueType =
typename GridType::ValueType;
3111 using Int32TreeType =
typename Int32GridType::TreeType;
3113 using BoolTreeType =
typename TreeType::template ValueConverter<bool>::Type;
3120 distGrid->setTransform(transform.
copy());
3122 ValueType exteriorWidth = ValueType(exteriorBandWidth);
3123 ValueType interiorWidth = ValueType(interiorBandWidth);
3127 if (!std::isfinite(exteriorWidth) || std::isnan(interiorWidth)) {
3128 std::stringstream msg;
3129 msg <<
"Illegal narrow band width: exterior = " << exteriorWidth
3130 <<
", interior = " << interiorWidth;
3135 const ValueType voxelSize = ValueType(transform.
voxelSize()[0]);
3137 if (!std::isfinite(voxelSize) ||
math::isZero(voxelSize)) {
3138 std::stringstream msg;
3139 msg <<
"Illegal transform, voxel size = " << voxelSize;
3145 exteriorWidth *= voxelSize;
3149 interiorWidth *= voxelSize;
3157 Int32GridType* indexGrid =
nullptr;
3159 typename Int32GridType::Ptr temporaryIndexGrid;
3161 if (polygonIndexGrid) {
3162 indexGrid = polygonIndexGrid;
3165 indexGrid = temporaryIndexGrid.get();
3168 indexGrid->newTree();
3169 indexGrid->setTransform(transform.
copy());
3171 if (computeSignedDistanceField) {
3175 interiorWidth = ValueType(0.0);
3178 TreeType& distTree = distGrid->tree();
3179 Int32TreeType& indexTree = indexGrid->tree();
3188 using DataTable = tbb::enumerable_thread_specific<typename VoxelizationDataType::Ptr>;
3194 const tbb::blocked_range<size_t> polygonRange(0, mesh.polygonCount());
3196 tbb::parallel_for(polygonRange, Voxelizer(data, mesh, &interrupter));
3198 for (
typename DataTable::iterator i = data.begin(); i != data.end(); ++i) {
3199 VoxelizationDataType& dataItem = **i;
3201 distTree, indexTree, dataItem.distTree, dataItem.indexTree);
3207 if (interrupter.wasInterrupted(30))
return distGrid;
3214 if (computeSignedDistanceField) {
3219 std::vector<LeafNodeType*> nodes;
3220 nodes.reserve(distTree.leafCount());
3221 distTree.getNodes(nodes);
3223 const tbb::blocked_range<size_t> nodeRange(0, nodes.size());
3228 tbb::parallel_for(nodeRange, SignOp(nodes, distTree, indexTree, mesh));
3230 if (interrupter.wasInterrupted(45))
return distGrid;
3233 if (removeIntersectingVoxels) {
3235 tbb::parallel_for(nodeRange,
3238 tbb::parallel_for(nodeRange,
3240 nodes, distTree, indexTree));
3247 if (interrupter.wasInterrupted(50))
return distGrid;
3249 if (distTree.activeVoxelCount() == 0) {
3251 distTree.root().setBackground(exteriorWidth,
false);
3257 std::vector<LeafNodeType*> nodes;
3258 nodes.reserve(distTree.leafCount());
3259 distTree.getNodes(nodes);
3261 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3263 nodes, voxelSize, !computeSignedDistanceField));
3267 if (computeSignedDistanceField) {
3268 distTree.root().setBackground(exteriorWidth,
false);
3274 if (interrupter.wasInterrupted(54))
return distGrid;
3281 const ValueType minBandWidth = voxelSize * ValueType(2.0);
3283 if (interiorWidth > minBandWidth || exteriorWidth > minBandWidth) {
3286 BoolTreeType maskTree(
false);
3289 std::vector<LeafNodeType*> nodes;
3290 nodes.reserve(distTree.leafCount());
3291 distTree.getNodes(nodes);
3294 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), op);
3300 float progress = 54.0f, step = 0.0f;
3302 2.0 * std::ceil((
std::max(interiorWidth, exteriorWidth) - minBandWidth) / voxelSize);
3304 if (estimated <
double(maxIterations)) {
3305 maxIterations = unsigned(estimated);
3306 step = 40.0f / float(maxIterations);
3309 std::vector<typename BoolTreeType::LeafNodeType*> maskNodes;
3314 if (interrupter.wasInterrupted(
int(progress)))
return distGrid;
3316 const size_t maskNodeCount = maskTree.leafCount();
3317 if (maskNodeCount == 0)
break;
3320 maskNodes.reserve(maskNodeCount);
3321 maskTree.getNodes(maskNodes);
3323 const tbb::blocked_range<size_t> range(0, maskNodes.size());
3325 tbb::parallel_for(range,
3329 mesh, exteriorWidth, interiorWidth, voxelSize);
3331 if ((++count) >= maxIterations)
break;
3336 if (interrupter.wasInterrupted(94))
return distGrid;
3338 if (!polygonIndexGrid) indexGrid->clear();
3346 if (computeSignedDistanceField && renormalizeValues) {
3348 std::vector<LeafNodeType*> nodes;
3349 nodes.reserve(distTree.leafCount());
3350 distTree.getNodes(nodes);
3352 std::unique_ptr<ValueType[]> buffer{
new ValueType[LeafNodeType::SIZE * nodes.size()]};
3354 const ValueType offset = ValueType(0.8 * voxelSize);
3356 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3359 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3361 distTree, nodes, buffer.get(), voxelSize));
3363 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3366 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3371 if (interrupter.wasInterrupted(99))
return distGrid;
3378 if (trimNarrowBand &&
std::min(interiorWidth, exteriorWidth) < voxelSize * ValueType(4.0)) {
3380 std::vector<LeafNodeType*> nodes;
3381 nodes.reserve(distTree.leafCount());
3382 distTree.getNodes(nodes);
3384 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3386 nodes, exteriorWidth, computeSignedDistanceField ? interiorWidth : exteriorWidth));
3389 distTree, exteriorWidth, computeSignedDistanceField ? -interiorWidth : -exteriorWidth);
3396 template <
typename Gr
idType,
typename MeshDataAdapter>
3397 inline typename GridType::Ptr
3401 float exteriorBandWidth,
3402 float interiorBandWidth,
3407 return meshToVolume<GridType>(nullInterrupter, mesh, transform,
3408 exteriorBandWidth, interiorBandWidth, flags, polygonIndexGrid);
3419 template<
typename Gr
idType,
typename Interrupter>
3420 inline typename std::enable_if<std::is_floating_point<typename GridType::ValueType>::value,
3421 typename GridType::Ptr>::type
3423 Interrupter& interrupter,
3424 const openvdb::math::Transform& xform,
3425 const std::vector<Vec3s>& points,
3426 const std::vector<Vec3I>& triangles,
3427 const std::vector<Vec4I>& quads,
3430 bool unsignedDistanceField =
false)
3432 if (points.empty()) {
3433 return typename GridType::Ptr(
new GridType(
typename GridType::ValueType(exBandWidth)));
3436 const size_t numPoints = points.size();
3437 std::unique_ptr<Vec3s[]> indexSpacePoints{
new Vec3s[numPoints]};
3440 tbb::parallel_for(tbb::blocked_range<size_t>(0, numPoints),
3441 mesh_to_volume_internal::TransformPoints<Vec3s>(
3442 &points[0], indexSpacePoints.get(), xform));
3446 if (quads.empty()) {
3448 QuadAndTriangleDataAdapter<Vec3s, Vec3I>
3449 mesh(indexSpacePoints.get(), numPoints, &triangles[0], triangles.size());
3451 return meshToVolume<GridType>(
3452 interrupter, mesh, xform, exBandWidth, inBandWidth, conversionFlags);
3454 }
else if (triangles.empty()) {
3456 QuadAndTriangleDataAdapter<Vec3s, Vec4I>
3457 mesh(indexSpacePoints.get(), numPoints, &quads[0], quads.size());
3459 return meshToVolume<GridType>(
3460 interrupter, mesh, xform, exBandWidth, inBandWidth, conversionFlags);
3465 const size_t numPrimitives = triangles.size() + quads.size();
3466 std::unique_ptr<Vec4I[]> prims{
new Vec4I[numPrimitives]};
3468 for (
size_t n = 0, N = triangles.size(); n < N; ++n) {
3469 const Vec3I& triangle = triangles[n];
3470 Vec4I& prim = prims[n];
3471 prim[0] = triangle[0];
3472 prim[1] = triangle[1];
3473 prim[2] = triangle[2];
3477 const size_t offset = triangles.size();
3478 for (
size_t n = 0, N = quads.size(); n < N; ++n) {
3479 prims[offset + n] = quads[n];
3482 QuadAndTriangleDataAdapter<Vec3s, Vec4I>
3483 mesh(indexSpacePoints.get(), numPoints, prims.get(), numPrimitives);
3485 return meshToVolume<GridType>(interrupter, mesh, xform,
3486 exBandWidth, inBandWidth, conversionFlags);
3492 template<
typename Gr
idType,
typename Interrupter>
3493 inline typename std::enable_if<!std::is_floating_point<typename GridType::ValueType>::value,
3494 typename GridType::Ptr>::type
3497 const math::Transform& ,
3498 const std::vector<Vec3s>& ,
3499 const std::vector<Vec3I>& ,
3500 const std::vector<Vec4I>& ,
3506 "mesh to volume conversion is supported only for scalar floating-point grids");
3516 template<
typename Gr
idType>
3517 inline typename GridType::Ptr
3519 const openvdb::math::Transform& xform,
3520 const std::vector<Vec3s>& points,
3521 const std::vector<Vec3I>& triangles,
3525 std::vector<Vec4I> quads(0);
3526 return doMeshConversion<GridType>(nullInterrupter, xform, points, triangles, quads,
3527 halfWidth, halfWidth);
3531 template<
typename Gr
idType,
typename Interrupter>
3532 inline typename GridType::Ptr
3534 Interrupter& interrupter,
3535 const openvdb::math::Transform& xform,
3536 const std::vector<Vec3s>& points,
3537 const std::vector<Vec3I>& triangles,
3540 std::vector<Vec4I> quads(0);
3541 return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3542 halfWidth, halfWidth);
3546 template<
typename Gr
idType>
3547 inline typename GridType::Ptr
3549 const openvdb::math::Transform& xform,
3550 const std::vector<Vec3s>& points,
3551 const std::vector<Vec4I>& quads,
3555 std::vector<Vec3I> triangles(0);
3556 return doMeshConversion<GridType>(nullInterrupter, xform, points, triangles, quads,
3557 halfWidth, halfWidth);
3561 template<
typename Gr
idType,
typename Interrupter>
3562 inline typename GridType::Ptr
3564 Interrupter& interrupter,
3565 const openvdb::math::Transform& xform,
3566 const std::vector<Vec3s>& points,
3567 const std::vector<Vec4I>& quads,
3570 std::vector<Vec3I> triangles(0);
3571 return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3572 halfWidth, halfWidth);
3576 template<
typename Gr
idType>
3577 inline typename GridType::Ptr
3579 const openvdb::math::Transform& xform,
3580 const std::vector<Vec3s>& points,
3581 const std::vector<Vec3I>& triangles,
3582 const std::vector<Vec4I>& quads,
3586 return doMeshConversion<GridType>(nullInterrupter, xform, points, triangles, quads,
3587 halfWidth, halfWidth);
3591 template<
typename Gr
idType,
typename Interrupter>
3592 inline typename GridType::Ptr
3594 Interrupter& interrupter,
3595 const openvdb::math::Transform& xform,
3596 const std::vector<Vec3s>& points,
3597 const std::vector<Vec3I>& triangles,
3598 const std::vector<Vec4I>& quads,
3601 return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3602 halfWidth, halfWidth);
3606 template<
typename Gr
idType>
3607 inline typename GridType::Ptr
3609 const openvdb::math::Transform& xform,
3610 const std::vector<Vec3s>& points,
3611 const std::vector<Vec3I>& triangles,
3612 const std::vector<Vec4I>& quads,
3617 return doMeshConversion<GridType>(nullInterrupter, xform, points, triangles,
3618 quads, exBandWidth, inBandWidth);
3622 template<
typename Gr
idType,
typename Interrupter>
3623 inline typename GridType::Ptr
3625 Interrupter& interrupter,
3626 const openvdb::math::Transform& xform,
3627 const std::vector<Vec3s>& points,
3628 const std::vector<Vec3I>& triangles,
3629 const std::vector<Vec4I>& quads,
3633 return doMeshConversion<GridType>(interrupter, xform, points, triangles,
3634 quads, exBandWidth, inBandWidth);
3638 template<
typename Gr
idType>
3639 inline typename GridType::Ptr
3641 const openvdb::math::Transform& xform,
3642 const std::vector<Vec3s>& points,
3643 const std::vector<Vec3I>& triangles,
3644 const std::vector<Vec4I>& quads,
3648 return doMeshConversion<GridType>(nullInterrupter, xform, points, triangles, quads,
3649 bandWidth, bandWidth,
true);
3653 template<
typename Gr
idType,
typename Interrupter>
3654 inline typename GridType::Ptr
3656 Interrupter& interrupter,
3657 const openvdb::math::Transform& xform,
3658 const std::vector<Vec3s>& points,
3659 const std::vector<Vec3I>& triangles,
3660 const std::vector<Vec4I>& quads,
3663 return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3664 bandWidth, bandWidth,
true);
3672 inline std::ostream&
3675 ostr <<
"{[ " << rhs.
mXPrim <<
", " << rhs.
mXDist <<
"]";
3676 ostr <<
" [ " << rhs.
mYPrim <<
", " << rhs.
mYDist <<
"]";
3677 ostr <<
" [ " << rhs.
mZPrim <<
", " << rhs.
mZDist <<
"]}";
3682 inline MeshToVoxelEdgeData::EdgeData
3697 const std::vector<Vec3s>& pointList,
3698 const std::vector<Vec4I>& polygonList);
3700 void run(
bool threaded =
true);
3703 inline void operator() (
const tbb::blocked_range<size_t> &range);
3711 struct Primitive {
Vec3d a, b, c, d;
Int32 index; };
3713 template<
bool IsQuad>
3714 inline void voxelize(
const Primitive&);
3716 template<
bool IsQuad>
3717 inline bool evalPrimitive(
const Coord&,
const Primitive&);
3719 inline bool rayTriangleIntersection(
const Vec3d& origin,
const Vec3d& dir,
3726 const std::vector<Vec3s>& mPointList;
3727 const std::vector<Vec4I>& mPolygonList;
3730 using IntTreeT = TreeType::ValueConverter<Int32>::Type;
3731 IntTreeT mLastPrimTree;
3737 MeshToVoxelEdgeData::GenEdgeData::GenEdgeData(
3738 const std::vector<Vec3s>& pointList,
3739 const std::vector<Vec4I>& polygonList)
3742 , mPointList(pointList)
3743 , mPolygonList(polygonList)
3745 , mLastPrimAccessor(mLastPrimTree)
3754 , mPointList(rhs.mPointList)
3755 , mPolygonList(rhs.mPolygonList)
3757 , mLastPrimAccessor(mLastPrimTree)
3766 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPolygonList.size()), *
this);
3768 (*this)(tbb::blocked_range<size_t>(0, mPolygonList.size()));
3777 using NodeChainType = RootNodeType::NodeChainType;
3778 static_assert(boost::mpl::size<NodeChainType>::value > 1,
"expected tree height > 1");
3779 using InternalNodeType = boost::mpl::at<NodeChainType, boost::mpl::int_<1> >::type;
3787 for ( ; leafIt; ++leafIt) {
3788 ijk = leafIt->origin();
3794 mAccessor.addLeaf(rhs.mAccessor.
probeLeaf(ijk));
3795 InternalNodeType* node = rhs.mAccessor.
getNode<InternalNodeType>();
3797 rhs.mAccessor.
clear();
3807 if (!lhsLeafPt->isValueOn(offset)) {
3808 lhsLeafPt->setValueOn(offset, rhsValue);
3811 EdgeData& lhsValue = const_cast<EdgeData&>(lhsLeafPt->getValue(offset));
3840 for (
size_t n = range.begin(); n < range.end(); ++n) {
3842 const Vec4I& verts = mPolygonList[n];
3844 prim.index =
Int32(n);
3845 prim.a =
Vec3d(mPointList[verts[0]]);
3846 prim.b =
Vec3d(mPointList[verts[1]]);
3847 prim.c =
Vec3d(mPointList[verts[2]]);
3850 prim.d =
Vec3d(mPointList[verts[3]]);
3851 voxelize<true>(prim);
3853 voxelize<false>(prim);
3859 template<
bool IsQuad>
3861 MeshToVoxelEdgeData::GenEdgeData::voxelize(
const Primitive& prim)
3863 std::deque<Coord> coordList;
3867 coordList.push_back(ijk);
3869 evalPrimitive<IsQuad>(ijk, prim);
3871 while (!coordList.empty()) {
3873 ijk = coordList.back();
3874 coordList.pop_back();
3876 for (
Int32 i = 0; i < 26; ++i) {
3879 if (prim.index != mLastPrimAccessor.getValue(nijk)) {
3880 mLastPrimAccessor.setValue(nijk, prim.index);
3881 if(evalPrimitive<IsQuad>(nijk, prim)) coordList.push_back(nijk);
3888 template<
bool IsQuad>
3890 MeshToVoxelEdgeData::GenEdgeData::evalPrimitive(
const Coord& ijk,
const Primitive& prim)
3892 Vec3d uvw, org(ijk[0], ijk[1], ijk[2]);
3893 bool intersecting =
false;
3897 mAccessor.probeValue(ijk, edgeData);
3900 double dist = (org -
3903 if (rayTriangleIntersection(org,
Vec3d(1.0, 0.0, 0.0), prim.a, prim.c, prim.b, t)) {
3904 if (t < edgeData.mXDist) {
3905 edgeData.mXDist = float(t);
3906 edgeData.mXPrim = prim.index;
3907 intersecting =
true;
3911 if (rayTriangleIntersection(org,
Vec3d(0.0, 1.0, 0.0), prim.a, prim.c, prim.b, t)) {
3912 if (t < edgeData.mYDist) {
3913 edgeData.mYDist = float(t);
3914 edgeData.mYPrim = prim.index;
3915 intersecting =
true;
3919 if (rayTriangleIntersection(org,
Vec3d(0.0, 0.0, 1.0), prim.a, prim.c, prim.b, t)) {
3920 if (t < edgeData.mZDist) {
3921 edgeData.mZDist = float(t);
3922 edgeData.mZPrim = prim.index;
3923 intersecting =
true;
3929 double secondDist = (org -
3932 if (secondDist < dist) dist = secondDist;
3934 if (rayTriangleIntersection(org,
Vec3d(1.0, 0.0, 0.0), prim.a, prim.d, prim.c, t)) {
3935 if (t < edgeData.mXDist) {
3936 edgeData.mXDist = float(t);
3937 edgeData.mXPrim = prim.index;
3938 intersecting =
true;
3942 if (rayTriangleIntersection(org,
Vec3d(0.0, 1.0, 0.0), prim.a, prim.d, prim.c, t)) {
3943 if (t < edgeData.mYDist) {
3944 edgeData.mYDist = float(t);
3945 edgeData.mYPrim = prim.index;
3946 intersecting =
true;
3950 if (rayTriangleIntersection(org,
Vec3d(0.0, 0.0, 1.0), prim.a, prim.d, prim.c, t)) {
3951 if (t < edgeData.mZDist) {
3952 edgeData.mZDist = float(t);
3953 edgeData.mZPrim = prim.index;
3954 intersecting =
true;
3959 if (intersecting) mAccessor.setValue(ijk, edgeData);
3961 return (dist < 0.86602540378443861);
3966 MeshToVoxelEdgeData::GenEdgeData::rayTriangleIntersection(
3977 double divisor = s1.
dot(e1);
3978 if (!(std::abs(divisor) > 0.0))
return false;
3982 double inv_divisor = 1.0 / divisor;
3983 Vec3d d = origin - a;
3984 double b1 = d.
dot(s1) * inv_divisor;
3986 if (b1 < 0.0 || b1 > 1.0)
return false;
3989 double b2 = dir.
dot(s2) * inv_divisor;
3991 if (b2 < 0.0 || (b1 + b2) > 1.0)
return false;
3995 t = e2.dot(s2) * inv_divisor;
3996 return (t < 0.0) ? false :
true;
4012 const std::vector<Vec3s>& pointList,
4013 const std::vector<Vec4I>& polygonList)
4027 std::vector<Vec3d>& points,
4028 std::vector<Index32>& primitives)
4038 point[0] = double(coord[0]) + data.
mXDist;
4039 point[1] = double(coord[1]);
4040 point[2] = double(coord[2]);
4042 points.push_back(point);
4043 primitives.push_back(data.
mXPrim);
4047 point[0] = double(coord[0]);
4048 point[1] = double(coord[1]) + data.
mYDist;
4049 point[2] = double(coord[2]);
4051 points.push_back(point);
4052 primitives.push_back(data.
mYPrim);
4056 point[0] = double(coord[0]);
4057 point[1] = double(coord[1]);
4058 point[2] = double(coord[2]) + data.
mZDist;
4060 points.push_back(point);
4061 primitives.push_back(data.
mZPrim);
4071 point[0] = double(coord[0]);
4072 point[1] = double(coord[1]) + data.
mYDist;
4073 point[2] = double(coord[2]);
4075 points.push_back(point);
4076 primitives.push_back(data.
mYPrim);
4080 point[0] = double(coord[0]);
4081 point[1] = double(coord[1]);
4082 point[2] = double(coord[2]) + data.
mZDist;
4084 points.push_back(point);
4085 primitives.push_back(data.
mZPrim);
4093 point[0] = double(coord[0]);
4094 point[1] = double(coord[1]) + data.
mYDist;
4095 point[2] = double(coord[2]);
4097 points.push_back(point);
4098 primitives.push_back(data.
mYPrim);
4107 point[0] = double(coord[0]) + data.
mXDist;
4108 point[1] = double(coord[1]);
4109 point[2] = double(coord[2]);
4111 points.push_back(point);
4112 primitives.push_back(data.
mXPrim);
4116 point[0] = double(coord[0]);
4117 point[1] = double(coord[1]) + data.
mYDist;
4118 point[2] = double(coord[2]);
4120 points.push_back(point);
4121 primitives.push_back(data.
mYPrim);
4131 point[0] = double(coord[0]) + data.
mXDist;
4132 point[1] = double(coord[1]);
4133 point[2] = double(coord[2]);
4135 points.push_back(point);
4136 primitives.push_back(data.
mXPrim);
4145 point[0] = double(coord[0]) + data.
mXDist;
4146 point[1] = double(coord[1]);
4147 point[2] = double(coord[2]);
4149 points.push_back(point);
4150 primitives.push_back(data.
mXPrim);
4154 point[0] = double(coord[0]);
4155 point[1] = double(coord[1]);
4156 point[2] = double(coord[2]) + data.
mZDist;
4158 points.push_back(point);
4159 primitives.push_back(data.
mZPrim);
4168 point[0] = double(coord[0]);
4169 point[1] = double(coord[1]);
4170 point[2] = double(coord[2]) + data.
mZDist;
4172 points.push_back(point);
4173 primitives.push_back(data.
mZPrim);
4179 template<
typename Gr
idType,
typename VecType>
4180 inline typename GridType::Ptr
4182 const openvdb::math::Transform& xform,
4183 typename VecType::ValueType halfWidth)
4189 points[0] =
Vec3s(pmin[0], pmin[1], pmin[2]);
4190 points[1] =
Vec3s(pmin[0], pmin[1], pmax[2]);
4191 points[2] =
Vec3s(pmax[0], pmin[1], pmax[2]);
4192 points[3] =
Vec3s(pmax[0], pmin[1], pmin[2]);
4193 points[4] =
Vec3s(pmin[0], pmax[1], pmin[2]);
4194 points[5] =
Vec3s(pmin[0], pmax[1], pmax[2]);
4195 points[6] =
Vec3s(pmax[0], pmax[1], pmax[2]);
4196 points[7] =
Vec3s(pmax[0], pmax[1], pmin[2]);
4199 faces[0] =
Vec4I(0, 1, 2, 3);
4200 faces[1] =
Vec4I(7, 6, 5, 4);
4201 faces[2] =
Vec4I(4, 5, 1, 0);
4202 faces[3] =
Vec4I(6, 7, 3, 2);
4203 faces[4] =
Vec4I(0, 3, 7, 4);
4204 faces[5] =
Vec4I(1, 5, 6, 2);
4208 return meshToVolume<GridType>(mesh, xform, halfWidth, halfWidth);
4216 #endif // OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED