43 #ifndef OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED 44 #define OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED 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> 65 #include <boost/integer_traits.hpp> 66 #include <boost/scoped_array.hpp> 74 #include <type_traits> 139 template <
typename Gr
idType,
typename MeshDataAdapter>
140 inline typename GridType::Ptr
144 float exteriorBandWidth = 3.0f,
145 float interiorBandWidth = 3.0f,
165 template <
typename Gr
idType,
typename MeshDataAdapter,
typename Interrupter>
166 inline typename GridType::Ptr
168 Interrupter& interrupter,
171 float exteriorBandWidth = 3.0f,
172 float interiorBandWidth = 3.0f,
190 template<
typename Po
intType,
typename PolygonType>
194 const std::vector<PolygonType>& polygons)
195 : mPointArray(points.empty() ? nullptr : &points[0])
196 , mPointArraySize(points.size())
197 , mPolygonArray(polygons.empty() ? nullptr : &polygons[0])
198 , mPolygonArraySize(polygons.size())
203 const PolygonType* polygonArray,
size_t polygonArraySize)
204 : mPointArray(pointArray)
205 , mPointArraySize(pointArraySize)
206 , mPolygonArray(polygonArray)
207 , mPolygonArraySize(polygonArraySize)
216 return (PolygonType::size == 3 || mPolygonArray[n][3] ==
util::INVALID_IDX) ? 3 : 4;
222 const PointType& p = mPointArray[mPolygonArray[n][int(v)]];
223 pos[0] = double(p[0]);
224 pos[1] = double(p[1]);
225 pos[2] = double(p[2]);
229 PointType
const *
const mPointArray;
230 size_t const mPointArraySize;
231 PolygonType
const *
const mPolygonArray;
232 size_t const mPolygonArraySize;
260 template<
typename Gr
idType>
261 inline typename GridType::Ptr
263 const openvdb::math::Transform& xform,
264 const std::vector<Vec3s>& points,
265 const std::vector<Vec3I>& triangles,
269 template<
typename Gr
idType,
typename Interrupter>
270 inline typename GridType::Ptr
272 Interrupter& interrupter,
273 const openvdb::math::Transform& xform,
274 const std::vector<Vec3s>& points,
275 const std::vector<Vec3I>& triangles,
294 template<
typename Gr
idType>
295 inline typename GridType::Ptr
297 const openvdb::math::Transform& xform,
298 const std::vector<Vec3s>& points,
299 const std::vector<Vec4I>& quads,
303 template<
typename Gr
idType,
typename Interrupter>
304 inline typename GridType::Ptr
306 Interrupter& interrupter,
307 const openvdb::math::Transform& xform,
308 const std::vector<Vec3s>& points,
309 const std::vector<Vec4I>& quads,
329 template<
typename Gr
idType>
330 inline typename GridType::Ptr
332 const openvdb::math::Transform& xform,
333 const std::vector<Vec3s>& points,
334 const std::vector<Vec3I>& triangles,
335 const std::vector<Vec4I>& quads,
339 template<
typename Gr
idType,
typename Interrupter>
340 inline typename GridType::Ptr
342 Interrupter& interrupter,
343 const openvdb::math::Transform& xform,
344 const std::vector<Vec3s>& points,
345 const std::vector<Vec3I>& triangles,
346 const std::vector<Vec4I>& quads,
368 template<
typename Gr
idType>
369 inline typename GridType::Ptr
371 const openvdb::math::Transform& xform,
372 const std::vector<Vec3s>& points,
373 const std::vector<Vec3I>& triangles,
374 const std::vector<Vec4I>& quads,
379 template<
typename Gr
idType,
typename Interrupter>
380 inline typename GridType::Ptr
382 Interrupter& interrupter,
383 const openvdb::math::Transform& xform,
384 const std::vector<Vec3s>& points,
385 const std::vector<Vec3I>& triangles,
386 const std::vector<Vec4I>& quads,
405 template<
typename Gr
idType>
406 inline typename GridType::Ptr
408 const openvdb::math::Transform& xform,
409 const std::vector<Vec3s>& points,
410 const std::vector<Vec3I>& triangles,
411 const std::vector<Vec4I>& quads,
415 template<
typename Gr
idType,
typename Interrupter>
416 inline typename GridType::Ptr
418 Interrupter& interrupter,
419 const openvdb::math::Transform& xform,
420 const std::vector<Vec3s>& points,
421 const std::vector<Vec3I>& triangles,
422 const std::vector<Vec4I>& quads,
435 template<
typename Gr
idType,
typename VecType>
436 inline typename GridType::Ptr
438 const openvdb::math::Transform& xform,
451 template <
typename FloatTreeT>
469 : mXDist(dist), mYDist(dist), mZDist(dist)
512 void convert(
const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList);
517 void getEdgeData(
Accessor& acc,
const Coord& ijk,
518 std::vector<Vec3d>& points, std::vector<Index32>& primitives);
537 namespace mesh_to_volume_internal {
539 template<
typename Po
intType>
544 : mPointsIn(pointsIn), mPointsOut(pointsOut), mXform(&xform)
548 void operator()(
const tbb::blocked_range<size_t>& range)
const {
552 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
554 const PointType& wsP = mPointsIn[n];
555 pos[0] = double(wsP[0]);
556 pos[1] = double(wsP[1]);
557 pos[2] = double(wsP[2]);
559 pos = mXform->worldToIndex(pos);
561 PointType& isP = mPointsOut[n];
562 isP[0] =
typename PointType::value_type(pos[0]);
563 isP[1] =
typename PointType::value_type(pos[1]);
564 isP[2] =
typename PointType::value_type(pos[2]);
574 template<
typename ValueType>
577 static ValueType
epsilon() {
return ValueType(1e-7); }
585 template<
typename TreeType>
597 : mDistTree(&lhsDistTree)
598 , mIdxTree(&lhsIdxTree)
599 , mRhsDistNodes(rhsDistNodes)
600 , mRhsIdxNodes(rhsIdxNodes)
604 void operator()(
const tbb::blocked_range<size_t>& range)
const {
609 using DistValueType =
typename LeafNodeType::ValueType;
610 using IndexValueType =
typename Int32LeafNodeType::ValueType;
612 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
614 const Coord& origin = mRhsDistNodes[n]->origin();
619 DistValueType* lhsDistData = lhsDistNode->buffer().data();
620 IndexValueType* lhsIdxData = lhsIdxNode->buffer().data();
622 const DistValueType* rhsDistData = mRhsDistNodes[n]->buffer().data();
623 const IndexValueType* rhsIdxData = mRhsIdxNodes[n]->buffer().data();
626 for (
Index32 offset = 0; offset < LeafNodeType::SIZE; ++offset) {
630 const DistValueType& lhsValue = lhsDistData[offset];
631 const DistValueType& rhsValue = rhsDistData[offset];
633 if (rhsValue < lhsValue) {
634 lhsDistNode->setValueOn(offset, rhsValue);
635 lhsIdxNode->setValueOn(offset, rhsIdxData[offset]);
637 lhsIdxNode->setValueOn(offset,
638 std::min(lhsIdxData[offset], rhsIdxData[offset]));
643 delete mRhsDistNodes[n];
644 delete mRhsIdxNodes[n];
650 TreeType *
const mDistTree;
661 template<
typename TreeType>
667 : mNodes(nodes.empty() ? nullptr : &nodes[0]), mCoordinates(coordinates)
671 void operator()(
const tbb::blocked_range<size_t>& range)
const {
672 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
673 Coord& origin =
const_cast<Coord&
>(mNodes[n]->origin());
674 mCoordinates[n] = origin;
675 origin[0] =
static_cast<int>(n);
684 template<
typename TreeType>
690 : mNodes(nodes.empty() ? nullptr : &nodes[0]), mCoordinates(coordinates)
694 void operator()(
const tbb::blocked_range<size_t>& range)
const {
695 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
696 Coord& origin =
const_cast<Coord&
>(mNodes[n]->origin());
697 origin[0] = mCoordinates[n][0];
706 template<
typename TreeType>
713 size_t* offsets,
size_t numNodes,
const CoordBBox& bbox)
715 , mCoordinates(coordinates)
717 , mNumNodes(numNodes)
727 void operator()(
const tbb::blocked_range<size_t>& range)
const {
729 size_t* offsetsNextX = mOffsets;
730 size_t* offsetsPrevX = mOffsets + mNumNodes;
731 size_t* offsetsNextY = mOffsets + mNumNodes * 2;
732 size_t* offsetsPrevY = mOffsets + mNumNodes * 3;
733 size_t* offsetsNextZ = mOffsets + mNumNodes * 4;
734 size_t* offsetsPrevZ = mOffsets + mNumNodes * 5;
739 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
740 const Coord& origin = mCoordinates[n];
741 offsetsNextX[n] = findNeighbourNode(acc, origin, Coord(LeafNodeType::DIM, 0, 0));
742 offsetsPrevX[n] = findNeighbourNode(acc, origin, Coord(-LeafNodeType::DIM, 0, 0));
743 offsetsNextY[n] = findNeighbourNode(acc, origin, Coord(0, LeafNodeType::DIM, 0));
744 offsetsPrevY[n] = findNeighbourNode(acc, origin, Coord(0, -LeafNodeType::DIM, 0));
745 offsetsNextZ[n] = findNeighbourNode(acc, origin, Coord(0, 0, LeafNodeType::DIM));
746 offsetsPrevZ[n] = findNeighbourNode(acc, origin, Coord(0, 0, -LeafNodeType::DIM));
751 const Coord& start,
const Coord& step)
const 753 Coord ijk = start + step;
754 CoordBBox bbox(mBBox);
756 while (bbox.isInside(ijk)) {
758 if (node)
return static_cast<size_t>(node->origin()[0]);
762 return boost::integer_traits<size_t>::const_max;
767 TreeType
const *
const mTree;
768 Coord
const *
const mCoordinates;
769 size_t *
const mOffsets;
771 const size_t mNumNodes;
772 const CoordBBox mBBox;
776 template<
typename TreeType>
779 enum { INVALID_OFFSET = boost::integer_traits<size_t>::const_max };
787 mLeafNodes.reserve(tree.leafCount());
788 tree.getNodes(mLeafNodes);
790 if (mLeafNodes.empty())
return;
793 tree.evalLeafBoundingBox(bbox);
795 const tbb::blocked_range<size_t> range(0, mLeafNodes.size());
799 boost::scoped_array<Coord> coordinates(
new Coord[mLeafNodes.size()]);
800 tbb::parallel_for(range,
804 mOffsets.reset(
new size_t[mLeafNodes.size() * 6]);
808 tree, coordinates.get(), mOffsets.get(), mLeafNodes.size(), bbox));
814 size_t size()
const {
return mLeafNodes.size(); }
816 std::vector<LeafNodeType*>&
nodes() {
return mLeafNodes; }
817 const std::vector<LeafNodeType*>&
nodes()
const {
return mLeafNodes; }
821 const size_t*
offsetsPrevX()
const {
return mOffsets.get() + mLeafNodes.size(); }
823 const size_t*
offsetsNextY()
const {
return mOffsets.get() + mLeafNodes.size() * 2; }
824 const size_t*
offsetsPrevY()
const {
return mOffsets.get() + mLeafNodes.size() * 3; }
826 const size_t*
offsetsNextZ()
const {
return mOffsets.get() + mLeafNodes.size() * 4; }
827 const size_t*
offsetsPrevZ()
const {
return mOffsets.get() + mLeafNodes.size() * 5; }
830 std::vector<LeafNodeType*> mLeafNodes;
831 boost::scoped_array<size_t> mOffsets;
835 template<
typename TreeType>
848 : mStartNodeIndices(startNodeIndices.empty() ? nullptr : &startNodeIndices[0])
849 , mConnectivity(&connectivity)
854 void operator()(
const tbb::blocked_range<size_t>& range)
const {
856 std::vector<LeafNodeType*>& nodes = mConnectivity->nodes();
859 size_t idxA = 0, idxB = 1;
862 const size_t* nextOffsets = mConnectivity->offsetsNextZ();
863 const size_t* prevOffsets = mConnectivity->offsetsPrevZ();
869 step = LeafNodeType::DIM;
871 nextOffsets = mConnectivity->offsetsNextY();
872 prevOffsets = mConnectivity->offsetsPrevY();
874 }
else if (mAxis ==
X_AXIS) {
878 step = LeafNodeType::DIM * LeafNodeType::DIM;
880 nextOffsets = mConnectivity->offsetsNextX();
881 prevOffsets = mConnectivity->offsetsPrevX();
889 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
891 size_t startOffset = mStartNodeIndices[n];
892 size_t lastOffset = startOffset;
896 for (a = 0; a < int(LeafNodeType::DIM); ++a) {
897 for (b = 0; b < int(LeafNodeType::DIM); ++b) {
899 pos = LeafNodeType::coordToOffset(ijk);
900 size_t offset = startOffset;
903 while ( offset != ConnectivityTable::INVALID_OFFSET &&
904 traceVoxelLine(*nodes[offset], pos, step) ) {
907 offset = nextOffsets[offset];
912 while (offset != ConnectivityTable::INVALID_OFFSET) {
914 offset = nextOffsets[offset];
919 pos += step * (LeafNodeType::DIM - 1);
920 while ( offset != ConnectivityTable::INVALID_OFFSET &&
921 traceVoxelLine(*nodes[offset], pos, -step)) {
922 offset = prevOffsets[offset];
934 bool isOutside =
true;
936 for (
Index i = 0; i < LeafNodeType::DIM; ++i) {
944 if (!(dist >
ValueType(0.75))) isOutside =
false;
957 size_t const *
const mStartNodeIndices;
964 template<
typename LeafNodeType>
968 using ValueType =
typename LeafNodeType::ValueType;
969 using Queue = std::deque<Index>;
972 ValueType* data = node.buffer().data();
976 for (
Index pos = 0; pos < LeafNodeType::SIZE; ++pos) {
977 if (data[pos] < 0.0) seedPoints.push_back(pos);
980 if (seedPoints.empty())
return;
983 for (Queue::iterator it = seedPoints.begin(); it != seedPoints.end(); ++it) {
984 ValueType& dist = data[*it];
991 Index pos(0), nextPos(0);
993 while (!seedPoints.empty()) {
995 pos = seedPoints.back();
996 seedPoints.pop_back();
998 ValueType& dist = data[pos];
1000 if (!(dist < ValueType(0.0))) {
1004 ijk = LeafNodeType::offsetToLocalCoord(pos);
1007 nextPos = pos - LeafNodeType::DIM * LeafNodeType::DIM;
1008 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1011 if (ijk[0] != (LeafNodeType::DIM - 1)) {
1012 nextPos = pos + LeafNodeType::DIM * LeafNodeType::DIM;
1013 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1017 nextPos = pos - LeafNodeType::DIM;
1018 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1021 if (ijk[1] != (LeafNodeType::DIM - 1)) {
1022 nextPos = pos + LeafNodeType::DIM;
1023 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1028 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1031 if (ijk[2] != (LeafNodeType::DIM - 1)) {
1033 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1040 template<
typename LeafNodeType>
1044 bool updatedNode =
false;
1046 using ValueType =
typename LeafNodeType::ValueType;
1047 ValueType* data = node.buffer().data();
1051 bool updatedSign =
true;
1052 while (updatedSign) {
1054 updatedSign =
false;
1056 for (
Index pos = 0; pos < LeafNodeType::SIZE; ++pos) {
1058 ValueType& dist = data[pos];
1060 if (!(dist < ValueType(0.0)) && dist > ValueType(0.75)) {
1062 ijk = LeafNodeType::offsetToLocalCoord(pos);
1065 if (ijk[2] != 0 && data[pos - 1] < ValueType(0.0)) {
1067 dist = ValueType(-dist);
1070 }
else if (ijk[2] != (LeafNodeType::DIM - 1) && data[pos + 1] < ValueType(0.0)) {
1072 dist = ValueType(-dist);
1075 }
else if (ijk[1] != 0 && data[pos - LeafNodeType::DIM] < ValueType(0.0)) {
1077 dist = ValueType(-dist);
1080 }
else if (ijk[1] != (LeafNodeType::DIM - 1)
1081 && data[pos + LeafNodeType::DIM] < ValueType(0.0))
1084 dist = ValueType(-dist);
1087 }
else if (ijk[0] != 0
1088 && data[pos - LeafNodeType::DIM * LeafNodeType::DIM] < ValueType(0.0))
1091 dist = ValueType(-dist);
1094 }
else if (ijk[0] != (LeafNodeType::DIM - 1)
1095 && data[pos + LeafNodeType::DIM * LeafNodeType::DIM] < ValueType(0.0))
1098 dist = ValueType(-dist);
1103 updatedNode |= updatedSign;
1110 template<
typename TreeType>
1118 : mNodes(nodes.empty() ? nullptr : &nodes[0])
1119 , mChangedNodeMask(changedNodeMask)
1124 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1125 if (mChangedNodeMask[n]) {
1127 mChangedNodeMask[n] =
scanFill(*mNodes[n]);
1137 template<
typename ValueType>
1140 FillArray(ValueType* array,
const ValueType v) : mArray(array), mValue(v) { }
1143 const ValueType v = mValue;
1144 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1154 template<
typename ValueType>
1156 fillArray(ValueType* array,
const ValueType val,
const size_t length)
1158 const auto grainSize = std::max<size_t>(
1159 length / tbb::task_scheduler_init::default_num_threads(), 1024);
1160 const tbb::blocked_range<size_t> range(0, length, grainSize);
1165 template<
typename TreeType>
1173 const bool* changedNodeMask,
bool* changedVoxelMask)
1174 : mNodes(nodes.empty() ? nullptr : &nodes[0])
1175 , mChangedNodeMask(changedNodeMask)
1176 , mChangedVoxelMask(changedVoxelMask)
1181 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1183 if (mChangedNodeMask[n]) {
1184 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1186 ValueType* data = mNodes[n]->buffer().data();
1188 for (
Index pos = 0; pos < LeafNodeType::SIZE; ++pos) {
1204 template<
typename TreeType>
1213 bool* changedNodeMask,
bool* nodeMask,
bool* changedVoxelMask)
1214 : mConnectivity(&connectivity)
1215 , mChangedNodeMask(changedNodeMask)
1216 , mNodeMask(nodeMask)
1217 , mChangedVoxelMask(changedVoxelMask)
1223 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1225 if (!mChangedNodeMask[n]) {
1227 bool changedValue =
false;
1229 changedValue |= processZ(n,
true);
1230 changedValue |= processZ(n,
false);
1232 changedValue |= processY(n,
true);
1233 changedValue |= processY(n,
false);
1235 changedValue |= processX(n,
true);
1236 changedValue |= processX(n,
false);
1238 mNodeMask[n] = changedValue;
1246 const size_t offset =
1247 firstFace ? mConnectivity->offsetsPrevZ()[n] : mConnectivity->offsetsNextZ()[n];
1248 if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1250 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1252 const ValueType* lhsData = mConnectivity->nodes()[n]->buffer().data();
1253 const ValueType* rhsData = mConnectivity->nodes()[offset]->buffer().data();
1255 const Index lastOffset = LeafNodeType::DIM - 1;
1256 const Index lhsOffset =
1257 firstFace ? 0 : lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1259 Index tmpPos(0), pos(0);
1260 bool changedValue =
false;
1262 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
1263 tmpPos = x << (2 * LeafNodeType::LOG2DIM);
1264 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
1265 pos = tmpPos + (y << LeafNodeType::LOG2DIM);
1267 if (lhsData[pos + lhsOffset] >
ValueType(0.75)) {
1268 if (rhsData[pos + rhsOffset] <
ValueType(0.0)) {
1269 changedValue =
true;
1270 mask[pos + lhsOffset] =
true;
1276 return changedValue;
1284 const size_t offset =
1285 firstFace ? mConnectivity->offsetsPrevY()[n] : mConnectivity->offsetsNextY()[n];
1286 if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1288 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1290 const ValueType* lhsData = mConnectivity->nodes()[n]->buffer().data();
1291 const ValueType* rhsData = mConnectivity->nodes()[offset]->buffer().data();
1293 const Index lastOffset = LeafNodeType::DIM * (LeafNodeType::DIM - 1);
1294 const Index lhsOffset =
1295 firstFace ? 0 : lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1297 Index tmpPos(0), pos(0);
1298 bool changedValue =
false;
1300 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
1301 tmpPos = x << (2 * LeafNodeType::LOG2DIM);
1302 for (
Index z = 0; z < LeafNodeType::DIM; ++z) {
1305 if (lhsData[pos + lhsOffset] >
ValueType(0.75)) {
1306 if (rhsData[pos + rhsOffset] <
ValueType(0.0)) {
1307 changedValue =
true;
1308 mask[pos + lhsOffset] =
true;
1314 return changedValue;
1322 const size_t offset =
1323 firstFace ? mConnectivity->offsetsPrevX()[n] : mConnectivity->offsetsNextX()[n];
1324 if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1326 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1328 const ValueType* lhsData = mConnectivity->nodes()[n]->buffer().data();
1329 const ValueType* rhsData = mConnectivity->nodes()[offset]->buffer().data();
1331 const Index lastOffset = LeafNodeType::DIM * LeafNodeType::DIM * (LeafNodeType::DIM-1);
1332 const Index lhsOffset =
1333 firstFace ? 0 : lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1335 Index tmpPos(0), pos(0);
1336 bool changedValue =
false;
1338 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
1339 tmpPos = y << LeafNodeType::LOG2DIM;
1340 for (
Index z = 0; z < LeafNodeType::DIM; ++z) {
1343 if (lhsData[pos + lhsOffset] >
ValueType(0.75)) {
1344 if (rhsData[pos + rhsOffset] <
ValueType(0.0)) {
1345 changedValue =
true;
1346 mask[pos + lhsOffset] =
true;
1352 return changedValue;
1367 template<
typename TreeType,
typename MeshDataAdapter>
1375 using LocalData = std::pair<boost::shared_array<Vec3d>, boost::shared_array<bool> >;
1379 std::vector<LeafNodeType*>& distNodes,
1380 const TreeType& distTree,
1383 : mDistNodes(distNodes.empty() ? nullptr : &distNodes[0])
1384 , mDistTree(&distTree)
1385 , mIndexTree(&indexTree)
1399 Index xPos(0), yPos(0);
1400 Coord ijk, nijk, nodeMin, nodeMax;
1401 Vec3d cp, xyz, nxyz, dir1, dir2;
1403 LocalData& localData = mLocalDataTable->local();
1405 boost::shared_array<Vec3d>& points = localData.first;
1406 if (!points) points.reset(
new Vec3d[LeafNodeType::SIZE * 2]);
1408 boost::shared_array<bool>& mask = localData.second;
1409 if (!mask) mask.reset(
new bool[LeafNodeType::SIZE]);
1412 typename LeafNodeType::ValueOnCIter it;
1414 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1420 const Int32* idxData = idxNode->buffer().data();
1422 nodeMin = node.origin();
1423 nodeMax = nodeMin.offsetBy(LeafNodeType::DIM - 1);
1426 memset(mask.get(), 0,
sizeof(bool) * LeafNodeType::SIZE);
1428 for (it = node.cbeginValueOn(); it; ++it) {
1429 Index pos = it.pos();
1432 if (dist < 0.0 || dist > 0.75)
continue;
1434 ijk = node.offsetToGlobalCoord(pos);
1436 xyz[0] = double(ijk[0]);
1437 xyz[1] = double(ijk[1]);
1438 xyz[2] = double(ijk[2]);
1444 bool flipSign =
false;
1446 for (nijk[0] = bbox.min()[0]; nijk[0] <= bbox.max()[0] && !flipSign; ++nijk[0]) {
1447 xPos = (nijk[0] & (LeafNodeType::DIM - 1u)) << (2 * LeafNodeType::LOG2DIM);
1448 for (nijk[1]=bbox.min()[1]; nijk[1] <= bbox.max()[1] && !flipSign; ++nijk[1]) {
1449 yPos = xPos + ((nijk[1] & (LeafNodeType::DIM-1u)) << LeafNodeType::LOG2DIM);
1450 for (nijk[2] = bbox.min()[2]; nijk[2] <= bbox.max()[2]; ++nijk[2]) {
1451 pos = yPos + (nijk[2] & (LeafNodeType::DIM - 1u));
1453 const Int32& polyIdx = idxData[pos];
1458 const Index pointIndex = pos * 2;
1464 nxyz[0] = double(nijk[0]);
1465 nxyz[1] = double(nijk[1]);
1466 nxyz[2] = double(nijk[2]);
1468 Vec3d& point = points[pointIndex];
1470 point = closestPoint(nxyz, polyIdx);
1472 Vec3d& direction = points[pointIndex + 1];
1473 direction = nxyz - point;
1477 dir1 = xyz - points[pointIndex];
1480 if (points[pointIndex + 1].dot(dir1) > 0.0) {
1491 for (
Int32 m = 0; m < 26; ++m) {
1494 if (!bbox.isInside(nijk) && distAcc.
probeValue(nijk, nval) && nval<-0.75) {
1495 nxyz[0] = double(nijk[0]);
1496 nxyz[1] = double(nijk[1]);
1497 nxyz[2] = double(nijk[2]);
1499 cp = closestPoint(nxyz, idxAcc.
getValue(nijk));
1507 if (dir2.dot(dir1) > 0.0) {
1523 Vec3d a, b, c, cp, uvw;
1525 const size_t polygon = size_t(polyIdx);
1526 mMesh->getIndexSpacePoint(polygon, 0, a);
1527 mMesh->getIndexSpacePoint(polygon, 1, b);
1528 mMesh->getIndexSpacePoint(polygon, 2, c);
1532 if (4 == mMesh->vertexCount(polygon)) {
1534 mMesh->getIndexSpacePoint(polygon, 3, b);
1538 if ((center - c).lengthSqr() < (center - cp).lengthSqr()) {
1548 TreeType
const *
const mDistTree;
1559 template<
typename LeafNodeType>
1563 using NodeT = LeafNodeType;
1565 const Coord ijk = NodeT::offsetToLocalCoord(pos);
1569 mask[0] = ijk[0] != (NodeT::DIM - 1);
1571 mask[1] = ijk[0] != 0;
1573 mask[2] = ijk[1] != (NodeT::DIM - 1);
1575 mask[3] = ijk[1] != 0;
1577 mask[4] = ijk[2] != (NodeT::DIM - 1);
1579 mask[5] = ijk[2] != 0;
1583 mask[6] = mask[0] && mask[5];
1585 mask[7] = mask[1] && mask[5];
1587 mask[8] = mask[0] && mask[4];
1589 mask[9] = mask[1] && mask[4];
1591 mask[10] = mask[0] && mask[2];
1593 mask[11] = mask[1] && mask[2];
1595 mask[12] = mask[0] && mask[3];
1597 mask[13] = mask[1] && mask[3];
1599 mask[14] = mask[3] && mask[4];
1601 mask[15] = mask[3] && mask[5];
1603 mask[16] = mask[2] && mask[4];
1605 mask[17] = mask[2] && mask[5];
1609 mask[18] = mask[1] && mask[3] && mask[5];
1611 mask[19] = mask[1] && mask[3] && mask[4];
1613 mask[20] = mask[0] && mask[3] && mask[4];
1615 mask[21] = mask[0] && mask[3] && mask[5];
1617 mask[22] = mask[1] && mask[2] && mask[5];
1619 mask[23] = mask[1] && mask[2] && mask[4];
1621 mask[24] = mask[0] && mask[2] && mask[4];
1623 mask[25] = mask[0] && mask[2] && mask[5];
1627 template<
typename Compare,
typename LeafNodeType>
1631 using NodeT = LeafNodeType;
1634 if (mask[5] && Compare::check(data[pos - 1]))
return true;
1636 if (mask[4] && Compare::check(data[pos + 1]))
return true;
1638 if (mask[3] && Compare::check(data[pos - NodeT::DIM]))
return true;
1640 if (mask[2] && Compare::check(data[pos + NodeT::DIM]))
return true;
1642 if (mask[1] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM]))
return true;
1644 if (mask[0] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM]))
return true;
1646 if (mask[6] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM]))
return true;
1648 if (mask[7] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - 1]))
return true;
1650 if (mask[8] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + 1]))
return true;
1652 if (mask[9] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + 1]))
return true;
1654 if (mask[10] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM]))
return true;
1656 if (mask[11] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM]))
return true;
1658 if (mask[12] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM]))
return true;
1660 if (mask[13] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM]))
return true;
1662 if (mask[14] && Compare::check(data[pos - NodeT::DIM + 1]))
return true;
1664 if (mask[15] && Compare::check(data[pos - NodeT::DIM - 1]))
return true;
1666 if (mask[16] && Compare::check(data[pos + NodeT::DIM + 1]))
return true;
1668 if (mask[17] && Compare::check(data[pos + NodeT::DIM - 1]))
return true;
1670 if (mask[18] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM - 1]))
return true;
1672 if (mask[19] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM + 1]))
return true;
1674 if (mask[20] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM + 1]))
return true;
1676 if (mask[21] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM - 1]))
return true;
1678 if (mask[22] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM - 1]))
return true;
1680 if (mask[23] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM + 1]))
return true;
1682 if (mask[24] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM + 1]))
return true;
1684 if (mask[25] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM - 1]))
return true;
1690 template<
typename Compare,
typename AccessorType>
1694 for (
Int32 m = 0; m < 26; ++m) {
1704 template<
typename TreeType>
1714 , mNodes(nodes.empty() ? nullptr : &nodes[0])
1721 bool neighbourMask[26];
1723 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1728 typename LeafNodeType::ValueOnCIter it;
1729 for (it = node.cbeginValueOn(); it; ++it) {
1731 const Index pos = it.pos();
1734 if (dist < 0.0 || dist > 0.75)
continue;
1737 maskNodeInternalNeighbours<LeafNodeType>(pos, neighbourMask);
1739 const bool hasNegativeNeighbour =
1740 checkNeighbours<IsNegative, LeafNodeType>(pos, data, neighbourMask) ||
1741 checkNeighbours<IsNegative>(node.offsetToGlobalCoord(pos), acc, neighbourMask);
1743 if (!hasNegativeNeighbour) {
1756 template<
typename TreeType>
1767 : mNodes(nodes.empty() ? nullptr : &nodes[0])
1768 , mDistTree(&distTree)
1769 , mIndexTree(&indexTree)
1777 bool neighbourMask[26];
1779 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1782 ValueType* data = distNode.buffer().data();
1784 typename Int32TreeType::LeafNodeType* idxNode =
1787 typename LeafNodeType::ValueOnCIter it;
1788 for (it = distNode.cbeginValueOn(); it; ++it) {
1790 const Index pos = it.pos();
1792 if (!(data[pos] > 0.75))
continue;
1795 maskNodeInternalNeighbours<LeafNodeType>(pos, neighbourMask);
1797 const bool hasBoundaryNeighbour =
1798 checkNeighbours<Comp, LeafNodeType>(pos, data, neighbourMask) ||
1799 checkNeighbours<Comp>(distNode.offsetToGlobalCoord(pos),distAcc,neighbourMask);
1801 if (!hasBoundaryNeighbour) {
1802 distNode.setValueOff(pos);
1803 idxNode->setValueOff(pos);
1818 template<
typename NodeType>
1825 using NodeMaskType =
typename NodeType::NodeMaskType;
1827 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1828 const_cast<NodeMaskType&
>(mNodes[n]->getChildMask()).setOff();
1836 template<
typename TreeType>
1840 using RootNodeType =
typename TreeType::RootNodeType;
1841 using NodeChainType =
typename RootNodeType::NodeChainType;
1842 using InternalNodeType =
typename boost::mpl::at<NodeChainType, boost::mpl::int_<1> >::type;
1844 std::vector<InternalNodeType*> nodes;
1845 tree.getNodes(nodes);
1847 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
1852 template<
typename TreeType>
1858 std::vector<LeafNodeType*>& overlappingNodes)
1859 : mLhsTree(&lhsTree)
1860 , mRhsTree(&rhsTree)
1861 , mNodes(&overlappingNodes)
1867 std::vector<LeafNodeType*> rhsLeafNodes;
1869 rhsLeafNodes.reserve(mRhsTree->leafCount());
1872 mRhsTree->stealNodes(rhsLeafNodes);
1876 for (
size_t n = 0, N = rhsLeafNodes.size(); n < N; ++n) {
1877 if (!acc.
probeLeaf(rhsLeafNodes[n]->origin())) {
1880 mNodes->push_back(rhsLeafNodes[n]);
1886 TreeType *
const mLhsTree;
1887 TreeType *
const mRhsTree;
1888 std::vector<LeafNodeType*> *
const mNodes;
1892 template<
typename DistTreeType,
typename IndexTreeType>
1895 DistTreeType& rhsDist, IndexTreeType& rhsIdx)
1897 using DistLeafNodeType =
typename DistTreeType::LeafNodeType;
1898 using IndexLeafNodeType =
typename IndexTreeType::LeafNodeType;
1900 std::vector<DistLeafNodeType*> overlappingDistNodes;
1901 std::vector<IndexLeafNodeType*> overlappingIdxNodes;
1904 tbb::task_group tasks;
1910 if (!overlappingDistNodes.empty() && !overlappingIdxNodes.empty()) {
1911 tbb::parallel_for(tbb::blocked_range<size_t>(0, overlappingDistNodes.size()),
1913 &overlappingDistNodes[0], &overlappingIdxNodes[0]));
1923 template<
typename TreeType>
1926 using Ptr = std::unique_ptr<VoxelizationData>;
1930 using UCharTreeType =
typename TreeType::template ValueConverter<unsigned char>::Type;
1941 , indexAcc(indexTree)
1942 , primIdTree(MaxPrimId)
1943 , primIdAcc(primIdTree)
1959 if (mPrimCount == MaxPrimId || primIdTree.leafCount() > 1000) {
1964 return mPrimCount++;
1969 enum { MaxPrimId = 100 };
1971 unsigned char mPrimCount;
1975 template<
typename TreeType,
typename MeshDataAdapter,
typename Interrupter = util::NullInterrupter>
1981 using DataTable = tbb::enumerable_thread_specific<typename VoxelizationDataType::Ptr>;
1985 Interrupter* interrupter =
nullptr)
1986 : mDataTable(&dataTable)
1988 , mInterrupter(interrupter)
1999 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2002 tbb::task::self().cancel_group_execution();
2006 const size_t numVerts = mMesh->vertexCount(n);
2009 if (numVerts == 3 || numVerts == 4) {
2011 prim.index =
Int32(n);
2013 mMesh->getIndexSpacePoint(n, 0, prim.a);
2014 mMesh->getIndexSpacePoint(n, 1, prim.b);
2015 mMesh->getIndexSpacePoint(n, 2, prim.c);
2017 evalTriangle(prim, *dataPtr);
2019 if (numVerts == 4) {
2020 mMesh->getIndexSpacePoint(n, 3, prim.b);
2021 evalTriangle(prim, *dataPtr);
2029 bool wasInterrupted()
const {
return mInterrupter && mInterrupter->wasInterrupted(); }
2031 struct Triangle {
Vec3d a, b, c;
Int32 index; };
2035 enum { POLYGON_LIMIT = 1000 };
2037 SubTask(
const Triangle& prim,
DataTable& dataTable,
2038 int subdivisionCount,
size_t polygonCount)
2039 : mLocalDataTable(&dataTable)
2041 , mSubdivisionCount(subdivisionCount)
2042 , mPolygonCount(polygonCount)
2046 void operator()()
const 2048 if (mSubdivisionCount <= 0 || mPolygonCount >= POLYGON_LIMIT) {
2053 voxelizeTriangle(mPrim, *dataPtr);
2056 spawnTasks(mPrim, *mLocalDataTable, mSubdivisionCount, mPolygonCount);
2061 Triangle
const mPrim;
2062 int const mSubdivisionCount;
2063 size_t const mPolygonCount;
2066 inline static int evalSubdivisionCount(
const Triangle& prim)
2068 const double ax = prim.a[0], bx = prim.b[0], cx = prim.c[0];
2071 const double ay = prim.a[1], by = prim.b[1], cy = prim.c[1];
2074 const double az = prim.a[2], bz = prim.b[2], cz = prim.c[2];
2077 return int(
std::max(dx,
std::max(dy, dz)) /
double(TreeType::LeafNodeType::DIM * 2));
2082 const size_t polygonCount = mMesh->polygonCount();
2083 const int subdivisionCount =
2084 polygonCount < SubTask::POLYGON_LIMIT ? evalSubdivisionCount(prim) : 0;
2086 if (subdivisionCount <= 0) {
2087 voxelizeTriangle(prim, data);
2089 spawnTasks(prim, *mDataTable, subdivisionCount, polygonCount);
2093 static void spawnTasks(
2094 const Triangle& mainPrim,
DataTable& dataTable,
int subdivisionCount,
size_t polygonCount)
2096 subdivisionCount -= 1;
2099 tbb::task_group tasks;
2101 const Vec3d ac = (mainPrim.a + mainPrim.c) * 0.5;
2102 const Vec3d bc = (mainPrim.b + mainPrim.c) * 0.5;
2103 const Vec3d ab = (mainPrim.a + mainPrim.b) * 0.5;
2106 prim.index = mainPrim.index;
2108 prim.a = mainPrim.a;
2111 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount));
2116 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount));
2119 prim.b = mainPrim.b;
2121 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount));
2125 prim.c = mainPrim.c;
2126 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount));
2133 std::deque<Coord> coordList;
2136 ijk = Coord::floor(prim.a);
2137 coordList.push_back(ijk);
2139 computeDistance(ijk, prim, data);
2144 while (!coordList.empty()) {
2145 ijk = coordList.back();
2146 coordList.pop_back();
2148 for (
Int32 i = 0; i < 26; ++i) {
2152 if(computeDistance(nijk, prim, data)) coordList.push_back(nijk);
2158 static bool computeDistance(
const Coord& ijk,
const Triangle& prim,
VoxelizationDataType& data)
2160 Vec3d uvw, voxelCenter(ijk[0], ijk[1], ijk[2]);
2162 using ValueType =
typename TreeType::ValueType;
2164 const ValueType dist = ValueType((voxelCenter -
2169 if (dist < oldDist) {
2178 return !(dist > 0.75);
2183 Interrupter *
const mInterrupter;
2190 template<
typename TreeType>
2196 using BoolTreeType =
typename TreeType::template ValueConverter<bool>::Type;
2200 std::vector<BoolLeafNodeType*>& lhsNodes)
2201 : mRhsTree(&rhsTree), mLhsNodes(lhsNodes.empty() ? nullptr : &lhsNodes[0])
2209 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2214 if (rhsNode) lhsNode->topologyDifference(*rhsNode,
false);
2219 TreeType
const *
const mRhsTree;
2224 template<
typename LeafNodeTypeA,
typename LeafNodeTypeB>
2228 : mNodesA(nodesA.empty() ? nullptr : &nodesA[0])
2229 , mNodesB(nodesB.empty() ? nullptr : &nodesB[0])
2234 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2235 mNodesA[n]->topologyUnion(*mNodesB[n]);
2240 LeafNodeTypeA **
const mNodesA;
2241 LeafNodeTypeB **
const mNodesB;
2245 template<
typename TreeType>
2250 using BoolTreeType =
typename TreeType::template ValueConverter<bool>::Type;
2254 std::vector<LeafNodeType*>& nodes)
2256 , mNodes(nodes.empty() ? nullptr : &nodes[0])
2257 , mLocalMaskTree(false)
2258 , mMaskTree(&maskTree)
2264 , mNodes(rhs.mNodes)
2265 , mLocalMaskTree(false)
2266 , mMaskTree(&mLocalMaskTree)
2272 using Iterator =
typename LeafNodeType::ValueOnCIter;
2277 Coord ijk, nijk, localCorod;
2280 for (
size_t n = range.begin(); n != range.end(); ++n) {
2284 CoordBBox bbox = node.getNodeBoundingBox();
2289 for (Iterator it = node.cbeginValueOn(); it; ++it) {
2290 ijk = it.getCoord();
2293 localCorod = LeafNodeType::offsetToLocalCoord(pos);
2295 if (localCorod[2] <
int(LeafNodeType::DIM - 1)) {
2297 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2299 nijk = ijk.offsetBy(0, 0, 1);
2303 if (localCorod[2] > 0) {
2305 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2307 nijk = ijk.offsetBy(0, 0, -1);
2311 if (localCorod[1] <
int(LeafNodeType::DIM - 1)) {
2312 npos = pos + LeafNodeType::DIM;
2313 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2315 nijk = ijk.offsetBy(0, 1, 0);
2319 if (localCorod[1] > 0) {
2320 npos = pos - LeafNodeType::DIM;
2321 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2323 nijk = ijk.offsetBy(0, -1, 0);
2327 if (localCorod[0] <
int(LeafNodeType::DIM - 1)) {
2328 npos = pos + LeafNodeType::DIM * LeafNodeType::DIM;
2329 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2331 nijk = ijk.offsetBy(1, 0, 0);
2335 if (localCorod[0] > 0) {
2336 npos = pos - LeafNodeType::DIM * LeafNodeType::DIM;
2337 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2339 nijk = ijk.offsetBy(-1, 0, 0);
2349 TreeType
const *
const mTree;
2358 template<
typename TreeType,
typename MeshDataAdapter>
2366 using BoolTreeType =
typename TreeType::template ValueConverter<bool>::Type;
2377 : idx(idx_), x(x_), y(y_), z(z_), dist(dist_)
2387 std::vector<BoolLeafNodeType*>& maskNodes,
2395 : mMaskNodes(maskNodes.empty() ? nullptr : &maskNodes[0])
2396 , mMaskTree(&maskTree)
2397 , mDistTree(&distTree)
2398 , mIndexTree(&indexTree)
2400 , mNewMaskTree(false)
2402 , mUpdatedDistNodes()
2404 , mUpdatedIndexNodes()
2405 , mExteriorBandWidth(exteriorBandWidth)
2406 , mInteriorBandWidth(interiorBandWidth)
2407 , mVoxelSize(voxelSize)
2412 : mMaskNodes(rhs.mMaskNodes)
2413 , mMaskTree(rhs.mMaskTree)
2414 , mDistTree(rhs.mDistTree)
2415 , mIndexTree(rhs.mIndexTree)
2417 , mNewMaskTree(false)
2419 , mUpdatedDistNodes()
2421 , mUpdatedIndexNodes()
2422 , mExteriorBandWidth(rhs.mExteriorBandWidth)
2423 , mInteriorBandWidth(rhs.mInteriorBandWidth)
2424 , mVoxelSize(rhs.mVoxelSize)
2430 mDistNodes.insert(mDistNodes.end(), rhs.mDistNodes.begin(), rhs.mDistNodes.end());
2431 mIndexNodes.insert(mIndexNodes.end(), rhs.mIndexNodes.begin(), rhs.mIndexNodes.end());
2433 mUpdatedDistNodes.insert(mUpdatedDistNodes.end(),
2434 rhs.mUpdatedDistNodes.begin(), rhs.mUpdatedDistNodes.end());
2436 mUpdatedIndexNodes.insert(mUpdatedIndexNodes.end(),
2437 rhs.mUpdatedIndexNodes.begin(), rhs.mUpdatedIndexNodes.end());
2439 mNewMaskTree.merge(rhs.mNewMaskTree);
2448 std::vector<Fragment> fragments;
2449 fragments.reserve(256);
2451 std::unique_ptr<LeafNodeType> newDistNodePt;
2452 std::unique_ptr<Int32LeafNodeType> newIndexNodePt;
2454 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2457 if (maskNode.isEmpty())
continue;
2461 const Coord& origin = maskNode.origin();
2466 assert(!distNodePt == !indexNodePt);
2468 bool usingNewNodes =
false;
2470 if (!distNodePt && !indexNodePt) {
2474 if (!newDistNodePt.get() && !newIndexNodePt.get()) {
2475 newDistNodePt.reset(
new LeafNodeType(origin, backgroundDist));
2479 if ((backgroundDist <
ValueType(0.0)) !=
2480 (newDistNodePt->getValue(0) <
ValueType(0.0))) {
2481 newDistNodePt->buffer().fill(backgroundDist);
2484 newDistNodePt->setOrigin(origin);
2485 newIndexNodePt->setOrigin(origin);
2488 distNodePt = newDistNodePt.get();
2489 indexNodePt = newIndexNodePt.get();
2491 usingNewNodes =
true;
2498 for (
typename BoolLeafNodeType::ValueOnIter it = maskNode.beginValueOn(); it; ++it) {
2499 bbox.expand(it.getCoord());
2504 gatherFragments(fragments, bbox, distAcc, indexAcc);
2509 bbox = maskNode.getNodeBoundingBox();
2511 bool updatedLeafNodes =
false;
2513 for (
typename BoolLeafNodeType::ValueOnIter it = maskNode.beginValueOn(); it; ++it) {
2515 const Coord ijk = it.getCoord();
2517 if (updateVoxel(ijk, 5, fragments, *distNodePt, *indexNodePt, &updatedLeafNodes)) {
2519 for (
Int32 i = 0; i < 6; ++i) {
2521 if (bbox.isInside(nijk)) {
2522 mask.setOn(BoolLeafNodeType::coordToOffset(nijk));
2528 for (
Int32 i = 6; i < 26; ++i) {
2530 if (bbox.isInside(nijk)) {
2531 mask.setOn(BoolLeafNodeType::coordToOffset(nijk));
2537 if (updatedLeafNodes) {
2540 mask -= indexNodePt->getValueMask();
2542 for (
typename NodeMaskType::OnIterator it = mask.beginOn(); it; ++it) {
2544 const Index pos = it.pos();
2545 const Coord ijk = maskNode.origin() + LeafNodeType::offsetToLocalCoord(pos);
2547 if (updateVoxel(ijk, 6, fragments, *distNodePt, *indexNodePt)) {
2548 for (
Int32 i = 0; i < 6; ++i) {
2555 if (usingNewNodes) {
2556 newDistNodePt->topologyUnion(*newIndexNodePt);
2557 mDistNodes.push_back(newDistNodePt.release());
2558 mIndexNodes.push_back(newIndexNodePt.release());
2560 mUpdatedDistNodes.push_back(distNodePt);
2561 mUpdatedIndexNodes.push_back(indexNodePt);
2581 gatherFragments(std::vector<Fragment>& fragments,
const CoordBBox& bbox,
2585 const Coord nodeMin = bbox.min() & ~(LeafNodeType::DIM - 1);
2586 const Coord nodeMax = bbox.max() & ~(LeafNodeType::DIM - 1);
2591 for (ijk[0] = nodeMin[0]; ijk[0] <= nodeMax[0]; ijk[0] += LeafNodeType::DIM) {
2592 for (ijk[1] = nodeMin[1]; ijk[1] <= nodeMax[1]; ijk[1] += LeafNodeType::DIM) {
2593 for (ijk[2] = nodeMin[2]; ijk[2] <= nodeMax[2]; ijk[2] += LeafNodeType::DIM) {
2597 ijk.offsetBy(LeafNodeType::DIM - 1));
2598 gatherFragments(fragments, region, *distleaf, *indexAcc.
probeLeaf(ijk));
2604 std::sort(fragments.begin(), fragments.end());
2608 gatherFragments(std::vector<Fragment>& fragments,
const CoordBBox& bbox,
2611 const typename LeafNodeType::NodeMaskType& mask = distLeaf.getValueMask();
2612 const ValueType* distData = distLeaf.buffer().data();
2613 const Int32* idxData = idxLeaf.buffer().data();
2615 for (
int x = bbox.min()[0]; x <= bbox.max()[0]; ++x) {
2616 const Index xPos = (x & (LeafNodeType::DIM - 1u)) << (2 * LeafNodeType::LOG2DIM);
2617 for (
int y = bbox.min()[1]; y <= bbox.max()[1]; ++y) {
2618 const Index yPos = xPos + ((y & (LeafNodeType::DIM - 1u)) << LeafNodeType::LOG2DIM);
2619 for (
int z = bbox.min()[2]; z <= bbox.max()[2]; ++z) {
2620 const Index pos = yPos + (z & (LeafNodeType::DIM - 1u));
2621 if (mask.isOn(pos)) {
2622 fragments.push_back(
Fragment(idxData[pos],x,y,z, std::abs(distData[pos])));
2632 computeDistance(
const Coord& ijk,
const Int32 manhattanLimit,
2633 const std::vector<Fragment>& fragments,
Int32& closestPrimIdx)
const 2635 Vec3d a, b, c, uvw, voxelCenter(ijk[0], ijk[1], ijk[2]);
2639 for (
size_t n = 0, N = fragments.size(); n < N; ++n) {
2641 const Fragment& fragment = fragments[n];
2642 if (lastIdx == fragment.
idx)
continue;
2644 const Int32 dx = std::abs(fragment.
x - ijk[0]);
2645 const Int32 dy = std::abs(fragment.
y - ijk[1]);
2646 const Int32 dz = std::abs(fragment.
z - ijk[2]);
2648 const Int32 manhattan = dx + dy + dz;
2649 if (manhattan > manhattanLimit)
continue;
2651 lastIdx = fragment.
idx;
2653 const size_t polygon = size_t(lastIdx);
2655 mMesh->getIndexSpacePoint(polygon, 0, a);
2656 mMesh->getIndexSpacePoint(polygon, 1, b);
2657 mMesh->getIndexSpacePoint(polygon, 2, c);
2659 primDist = (voxelCenter -
2663 if (4 == mMesh->vertexCount(polygon)) {
2665 mMesh->getIndexSpacePoint(polygon, 3, b);
2668 a, b, c, voxelCenter, uvw)).lengthSqr();
2670 if (tmpDist < primDist) primDist = tmpDist;
2673 if (primDist < dist) {
2675 closestPrimIdx = lastIdx;
2679 return ValueType(std::sqrt(dist)) * mVoxelSize;
2685 updateVoxel(
const Coord& ijk,
const Int32 manhattanLimit,
2686 const std::vector<Fragment>& fragments,
2689 Int32 closestPrimIdx = 0;
2690 const ValueType distance = computeDistance(ijk, manhattanLimit, fragments, closestPrimIdx);
2692 const Index pos = LeafNodeType::coordToOffset(ijk);
2693 const bool inside = distLeaf.getValue(pos) <
ValueType(0.0);
2695 bool activateNeighbourVoxels =
false;
2697 if (!inside && distance < mExteriorBandWidth) {
2698 if (updatedLeafNodes) *updatedLeafNodes =
true;
2699 activateNeighbourVoxels = (distance + mVoxelSize) < mExteriorBandWidth;
2700 distLeaf.setValueOnly(pos, distance);
2701 idxLeaf.setValueOn(pos, closestPrimIdx);
2702 }
else if (inside && distance < mInteriorBandWidth) {
2703 if (updatedLeafNodes) *updatedLeafNodes =
true;
2704 activateNeighbourVoxels = (distance + mVoxelSize) < mInteriorBandWidth;
2705 distLeaf.setValueOnly(pos, -distance);
2706 idxLeaf.setValueOn(pos, closestPrimIdx);
2709 return activateNeighbourVoxels;
2716 TreeType *
const mDistTree;
2723 std::vector<LeafNodeType*> mDistNodes, mUpdatedDistNodes;
2724 std::vector<Int32LeafNodeType*> mIndexNodes, mUpdatedIndexNodes;
2726 const ValueType mExteriorBandWidth, mInteriorBandWidth, mVoxelSize;
2730 template<
typename TreeType>
2734 AddNodes(TreeType& tree, std::vector<LeafNodeType*>& nodes)
2735 : mTree(&tree) , mNodes(&nodes)
2741 std::vector<LeafNodeType*>& nodes = *mNodes;
2742 for (
size_t n = 0, N = nodes.size(); n < N; ++n) {
2752 template<
typename TreeType,
typename Int32TreeType,
typename BoolTreeType,
typename MeshDataAdapter>
2756 Int32TreeType& indexTree,
2757 BoolTreeType& maskTree,
2758 std::vector<typename BoolTreeType::LeafNodeType*>& maskNodes,
2760 typename TreeType::ValueType exteriorBandWidth,
2761 typename TreeType::ValueType interiorBandWidth,
2762 typename TreeType::ValueType voxelSize)
2765 distTree, indexTree, mesh, exteriorBandWidth, interiorBandWidth, voxelSize);
2767 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, maskNodes.size()), expandOp);
2769 tbb::parallel_for(tbb::blocked_range<size_t>(0, expandOp.
updatedIndexNodes().size()),
2773 tbb::task_group tasks;
2787 template<
typename TreeType>
2796 , mVoxelSize(voxelSize)
2797 , mUnsigned(unsignedDist)
2803 typename LeafNodeType::ValueOnIter iter;
2805 const bool udf = mUnsigned;
2806 const ValueType w[2] = { -mVoxelSize, mVoxelSize };
2808 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2810 for (iter = mNodes[n]->beginValueOn(); iter; ++iter) {
2812 val = w[udf || (val <
ValueType(0.0))] * std::sqrt(std::abs(val));
2820 const bool mUnsigned;
2825 template<
typename TreeType>
2833 : mNodes(nodes.empty() ? nullptr : &nodes[0])
2834 , mExBandWidth(exBandWidth)
2835 , mInBandWidth(inBandWidth)
2841 typename LeafNodeType::ValueOnIter iter;
2845 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2847 for (iter = mNodes[n]->beginValueOn(); iter; ++iter) {
2851 const bool inside = val <
ValueType(0.0);
2853 if (inside && !(val > inVal)) {
2856 }
else if (!inside && !(val < exVal)) {
2866 const ValueType mExBandWidth, mInBandWidth;
2870 template<
typename TreeType>
2877 : mNodes(nodes.empty() ? nullptr : &nodes[0]), mOffset(offset)
2885 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2887 typename LeafNodeType::ValueOnIter iter = mNodes[n]->beginValueOn();
2889 for (; iter; ++iter) {
2902 template<
typename TreeType>
2908 Renormalize(
const TreeType& tree,
const std::vector<LeafNodeType*>& nodes,
2911 , mNodes(nodes.empty() ? nullptr : &nodes[0])
2913 , mVoxelSize(voxelSize)
2928 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2930 ValueType* bufferData = &mBuffer[n * LeafNodeType::SIZE];
2932 typename LeafNodeType::ValueOnCIter iter = mNodes[n]->cbeginValueOn();
2933 for (; iter; ++iter) {
2937 ijk = iter.getCoord();
2939 up[0] = acc.
getValue(ijk.offsetBy(1, 0, 0)) - phi0;
2940 up[1] = acc.
getValue(ijk.offsetBy(0, 1, 0)) - phi0;
2941 up[2] = acc.
getValue(ijk.offsetBy(0, 0, 1)) - phi0;
2943 down[0] = phi0 - acc.
getValue(ijk.offsetBy(-1, 0, 0));
2944 down[1] = phi0 - acc.
getValue(ijk.offsetBy(0, -1, 0));
2945 down[2] = phi0 - acc.
getValue(ijk.offsetBy(0, 0, -1));
2952 bufferData[iter.pos()] = phi0 - dx * S * diff;
2958 TreeType
const *
const mTree;
2966 template<
typename TreeType>
2973 : mNodes(nodes.empty() ? nullptr : &nodes[0]), mBuffer(buffer)
2979 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2981 const ValueType* bufferData = &mBuffer[n * LeafNodeType::SIZE];
2983 typename LeafNodeType::ValueOnIter iter = mNodes[n]->beginValueOn();
2985 for (; iter; ++iter) {
2987 val =
std::min(val, bufferData[iter.pos()]);
3006 template <
typename FloatTreeT>
3012 ConnectivityTable nodeConnectivity(tree);
3014 std::vector<size_t> zStartNodes, yStartNodes, xStartNodes;
3016 for (
size_t n = 0; n < nodeConnectivity.size(); ++n) {
3017 if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevX()[n]) {
3018 xStartNodes.push_back(n);
3021 if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevY()[n]) {
3022 yStartNodes.push_back(n);
3025 if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevZ()[n]) {
3026 zStartNodes.push_back(n);
3032 tbb::parallel_for(tbb::blocked_range<size_t>(0, zStartNodes.size()),
3035 tbb::parallel_for(tbb::blocked_range<size_t>(0, yStartNodes.size()),
3038 tbb::parallel_for(tbb::blocked_range<size_t>(0, xStartNodes.size()),
3041 const size_t numLeafNodes = nodeConnectivity.size();
3042 const size_t numVoxels = numLeafNodes * FloatTreeT::LeafNodeType::SIZE;
3044 boost::scoped_array<bool> changedNodeMaskA(
new bool[numLeafNodes]);
3045 boost::scoped_array<bool> changedNodeMaskB(
new bool[numLeafNodes]);
3046 boost::scoped_array<bool> changedVoxelMask(
new bool[numVoxels]);
3052 const tbb::blocked_range<size_t> nodeRange(0, numLeafNodes);
3054 bool nodesUpdated =
false;
3057 nodeConnectivity.nodes(), changedNodeMaskA.get()));
3060 nodeConnectivity, changedNodeMaskA.get(), changedNodeMaskB.get(),
3061 changedVoxelMask.get()));
3063 changedNodeMaskA.swap(changedNodeMaskB);
3065 nodesUpdated =
false;
3066 for (
size_t n = 0; n < numLeafNodes; ++n) {
3067 nodesUpdated |= changedNodeMaskA[n];
3068 if (nodesUpdated)
break;
3073 nodeConnectivity.nodes(), changedNodeMaskA.get(), changedVoxelMask.get()));
3075 }
while (nodesUpdated);
3083 template <
typename Gr
idType,
typename MeshDataAdapter,
typename Interrupter>
3084 inline typename GridType::Ptr
3086 Interrupter& interrupter,
3089 float exteriorBandWidth,
3090 float interiorBandWidth,
3094 using GridTypePtr =
typename GridType::Ptr;
3095 using TreeType =
typename GridType::TreeType;
3096 using LeafNodeType =
typename TreeType::LeafNodeType;
3097 using ValueType =
typename GridType::ValueType;
3100 using Int32TreeType =
typename Int32GridType::TreeType;
3102 using BoolTreeType =
typename TreeType::template ValueConverter<bool>::Type;
3109 distGrid->setTransform(transform.
copy());
3111 ValueType exteriorWidth = ValueType(exteriorBandWidth);
3112 ValueType interiorWidth = ValueType(interiorBandWidth);
3116 if (!std::isfinite(exteriorWidth) || std::isnan(interiorWidth)) {
3117 std::stringstream msg;
3118 msg <<
"Illegal narrow band width: exterior = " << exteriorWidth
3119 <<
", interior = " << interiorWidth;
3124 const ValueType voxelSize = ValueType(transform.
voxelSize()[0]);
3126 if (!std::isfinite(voxelSize) ||
math::isZero(voxelSize)) {
3127 std::stringstream msg;
3128 msg <<
"Illegal transform, voxel size = " << voxelSize;
3134 exteriorWidth *= voxelSize;
3138 interiorWidth *= voxelSize;
3146 Int32GridType* indexGrid =
nullptr;
3148 typename Int32GridType::Ptr temporaryIndexGrid;
3150 if (polygonIndexGrid) {
3151 indexGrid = polygonIndexGrid;
3154 indexGrid = temporaryIndexGrid.get();
3157 indexGrid->newTree();
3158 indexGrid->setTransform(transform.
copy());
3160 if (computeSignedDistanceField) {
3164 interiorWidth = ValueType(0.0);
3167 TreeType& distTree = distGrid->tree();
3168 Int32TreeType& indexTree = indexGrid->tree();
3177 using DataTable = tbb::enumerable_thread_specific<typename VoxelizationDataType::Ptr>;
3183 const tbb::blocked_range<size_t> polygonRange(0, mesh.polygonCount());
3185 tbb::parallel_for(polygonRange, Voxelizer(data, mesh, &interrupter));
3187 for (
typename DataTable::iterator i = data.begin(); i != data.end(); ++i) {
3188 VoxelizationDataType& dataItem = **i;
3190 distTree, indexTree, dataItem.distTree, dataItem.indexTree);
3196 if (interrupter.wasInterrupted(30))
return distGrid;
3203 if (computeSignedDistanceField) {
3208 std::vector<LeafNodeType*> nodes;
3209 nodes.reserve(distTree.leafCount());
3210 distTree.getNodes(nodes);
3212 const tbb::blocked_range<size_t> nodeRange(0, nodes.size());
3217 tbb::parallel_for(nodeRange, SignOp(nodes, distTree, indexTree, mesh));
3219 if (interrupter.wasInterrupted(45))
return distGrid;
3222 if (removeIntersectingVoxels) {
3224 tbb::parallel_for(nodeRange,
3227 tbb::parallel_for(nodeRange,
3229 nodes, distTree, indexTree));
3236 if (interrupter.wasInterrupted(50))
return distGrid;
3238 if (distTree.activeVoxelCount() == 0) {
3240 distTree.root().setBackground(exteriorWidth,
false);
3246 std::vector<LeafNodeType*> nodes;
3247 nodes.reserve(distTree.leafCount());
3248 distTree.getNodes(nodes);
3250 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3252 nodes, voxelSize, !computeSignedDistanceField));
3256 if (computeSignedDistanceField) {
3257 distTree.root().setBackground(exteriorWidth,
false);
3263 if (interrupter.wasInterrupted(54))
return distGrid;
3270 const ValueType minBandWidth = voxelSize * ValueType(2.0);
3272 if (interiorWidth > minBandWidth || exteriorWidth > minBandWidth) {
3275 BoolTreeType maskTree(
false);
3278 std::vector<LeafNodeType*> nodes;
3279 nodes.reserve(distTree.leafCount());
3280 distTree.getNodes(nodes);
3283 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), op);
3289 float progress = 54.0f, step = 0.0f;
3291 2.0 * std::ceil((
std::max(interiorWidth, exteriorWidth) - minBandWidth) / voxelSize);
3293 if (estimated <
double(maxIterations)) {
3294 maxIterations = unsigned(estimated);
3295 step = 40.0f / float(maxIterations);
3298 std::vector<typename BoolTreeType::LeafNodeType*> maskNodes;
3303 if (interrupter.wasInterrupted(
int(progress)))
return distGrid;
3305 const size_t maskNodeCount = maskTree.leafCount();
3306 if (maskNodeCount == 0)
break;
3309 maskNodes.reserve(maskNodeCount);
3310 maskTree.getNodes(maskNodes);
3312 const tbb::blocked_range<size_t> range(0, maskNodes.size());
3314 tbb::parallel_for(range,
3318 mesh, exteriorWidth, interiorWidth, voxelSize);
3320 if ((++count) >= maxIterations)
break;
3325 if (interrupter.wasInterrupted(94))
return distGrid;
3327 if (!polygonIndexGrid) indexGrid->clear();
3335 if (computeSignedDistanceField && renormalizeValues) {
3337 std::vector<LeafNodeType*> nodes;
3338 nodes.reserve(distTree.leafCount());
3339 distTree.getNodes(nodes);
3341 boost::scoped_array<ValueType> buffer(
new ValueType[LeafNodeType::SIZE * nodes.size()]);
3343 const ValueType offset = ValueType(0.8 * voxelSize);
3345 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3348 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3350 distTree, nodes, buffer.get(), voxelSize));
3352 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3355 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3360 if (interrupter.wasInterrupted(99))
return distGrid;
3367 if (trimNarrowBand &&
std::min(interiorWidth, exteriorWidth) < voxelSize * ValueType(4.0)) {
3369 std::vector<LeafNodeType*> nodes;
3370 nodes.reserve(distTree.leafCount());
3371 distTree.getNodes(nodes);
3373 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3375 nodes, exteriorWidth, computeSignedDistanceField ? interiorWidth : exteriorWidth));
3378 distTree, exteriorWidth, computeSignedDistanceField ? -interiorWidth : -exteriorWidth);
3385 template <
typename Gr
idType,
typename MeshDataAdapter>
3386 inline typename GridType::Ptr
3390 float exteriorBandWidth,
3391 float interiorBandWidth,
3396 return meshToVolume<GridType>(nullInterrupter, mesh, transform,
3397 exteriorBandWidth, interiorBandWidth, flags, polygonIndexGrid);
3405 template<
typename Gr
idType,
typename Interrupter>
3406 inline typename std::enable_if<std::is_floating_point<typename GridType::ValueType>::value,
3407 typename GridType::Ptr>::type
3409 Interrupter& interrupter,
3410 const openvdb::math::Transform& xform,
3411 const std::vector<Vec3s>& points,
3412 const std::vector<Vec3I>& triangles,
3413 const std::vector<Vec4I>& quads,
3416 bool unsignedDistanceField =
false)
3418 if (points.empty()) {
3419 return typename GridType::Ptr(
new GridType(
typename GridType::ValueType(exBandWidth)));
3422 const size_t numPoints = points.size();
3423 boost::scoped_array<Vec3s> indexSpacePoints(
new Vec3s[numPoints]);
3426 tbb::parallel_for(tbb::blocked_range<size_t>(0, numPoints),
3428 &points[0], indexSpacePoints.get(), xform));
3432 if (quads.empty()) {
3435 mesh(indexSpacePoints.get(), numPoints, &triangles[0], triangles.size());
3437 return meshToVolume<GridType>(mesh, xform, exBandWidth, inBandWidth, conversionFlags);
3439 }
else if (triangles.empty()) {
3442 mesh(indexSpacePoints.get(), numPoints, &quads[0], quads.size());
3444 return meshToVolume<GridType>(mesh, xform, exBandWidth, inBandWidth, conversionFlags);
3449 const size_t numPrimitives = triangles.size() + quads.size();
3450 boost::scoped_array<Vec4I> prims(
new Vec4I[numPrimitives]);
3452 for (
size_t n = 0, N = triangles.size(); n < N; ++n) {
3453 const Vec3I& triangle = triangles[n];
3454 Vec4I& prim = prims[n];
3455 prim[0] = triangle[0];
3456 prim[1] = triangle[1];
3457 prim[2] = triangle[2];
3461 const size_t offset = triangles.size();
3462 for (
size_t n = 0, N = quads.size(); n < N; ++n) {
3463 prims[offset + n] = quads[n];
3467 mesh(indexSpacePoints.get(), numPoints, prims.get(), numPrimitives);
3469 return meshToVolume<GridType>(interrupter, mesh, xform,
3470 exBandWidth, inBandWidth, conversionFlags);
3476 template<
typename Gr
idType,
typename Interrupter>
3477 inline typename std::enable_if<!std::is_floating_point<typename GridType::ValueType>::value,
3478 typename GridType::Ptr>::type
3482 const std::vector<Vec3s>& ,
3483 const std::vector<Vec3I>& ,
3484 const std::vector<Vec4I>& ,
3490 "mesh to volume conversion is supported only for scalar floating-point grids");
3497 template<
typename Gr
idType>
3498 inline typename GridType::Ptr
3500 const openvdb::math::Transform& xform,
3501 const std::vector<Vec3s>& points,
3502 const std::vector<Vec3I>& triangles,
3506 std::vector<Vec4I> quads(0);
3507 return doMeshConversion<GridType>(nullInterrupter, xform, points, triangles, quads,
3508 halfWidth, halfWidth);
3512 template<
typename Gr
idType,
typename Interrupter>
3513 inline typename GridType::Ptr
3515 Interrupter& interrupter,
3516 const openvdb::math::Transform& xform,
3517 const std::vector<Vec3s>& points,
3518 const std::vector<Vec3I>& triangles,
3521 std::vector<Vec4I> quads(0);
3522 return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3523 halfWidth, halfWidth);
3527 template<
typename Gr
idType>
3528 inline typename GridType::Ptr
3530 const openvdb::math::Transform& xform,
3531 const std::vector<Vec3s>& points,
3532 const std::vector<Vec4I>& quads,
3536 std::vector<Vec3I> triangles(0);
3537 return doMeshConversion<GridType>(nullInterrupter, xform, points, triangles, quads,
3538 halfWidth, halfWidth);
3542 template<
typename Gr
idType,
typename Interrupter>
3543 inline typename GridType::Ptr
3545 Interrupter& interrupter,
3546 const openvdb::math::Transform& xform,
3547 const std::vector<Vec3s>& points,
3548 const std::vector<Vec4I>& quads,
3551 std::vector<Vec3I> triangles(0);
3552 return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3553 halfWidth, halfWidth);
3557 template<
typename Gr
idType>
3558 inline typename GridType::Ptr
3560 const openvdb::math::Transform& xform,
3561 const std::vector<Vec3s>& points,
3562 const std::vector<Vec3I>& triangles,
3563 const std::vector<Vec4I>& quads,
3567 return doMeshConversion<GridType>(nullInterrupter, xform, points, triangles, quads,
3568 halfWidth, halfWidth);
3572 template<
typename Gr
idType,
typename Interrupter>
3573 inline typename GridType::Ptr
3575 Interrupter& interrupter,
3576 const openvdb::math::Transform& xform,
3577 const std::vector<Vec3s>& points,
3578 const std::vector<Vec3I>& triangles,
3579 const std::vector<Vec4I>& quads,
3582 return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3583 halfWidth, halfWidth);
3587 template<
typename Gr
idType>
3588 inline typename GridType::Ptr
3590 const openvdb::math::Transform& xform,
3591 const std::vector<Vec3s>& points,
3592 const std::vector<Vec3I>& triangles,
3593 const std::vector<Vec4I>& quads,
3598 return doMeshConversion<GridType>(nullInterrupter, xform, points, triangles,
3599 quads, exBandWidth, inBandWidth);
3603 template<
typename Gr
idType,
typename Interrupter>
3604 inline typename GridType::Ptr
3606 Interrupter& interrupter,
3607 const openvdb::math::Transform& xform,
3608 const std::vector<Vec3s>& points,
3609 const std::vector<Vec3I>& triangles,
3610 const std::vector<Vec4I>& quads,
3614 return doMeshConversion<GridType>(interrupter, xform, points, triangles,
3615 quads, exBandWidth, inBandWidth);
3619 template<
typename Gr
idType>
3620 inline typename GridType::Ptr
3622 const openvdb::math::Transform& xform,
3623 const std::vector<Vec3s>& points,
3624 const std::vector<Vec3I>& triangles,
3625 const std::vector<Vec4I>& quads,
3629 return doMeshConversion<GridType>(nullInterrupter, xform, points, triangles, quads,
3630 bandWidth, bandWidth,
true);
3634 template<
typename Gr
idType,
typename Interrupter>
3635 inline typename GridType::Ptr
3637 Interrupter& interrupter,
3638 const openvdb::math::Transform& xform,
3639 const std::vector<Vec3s>& points,
3640 const std::vector<Vec3I>& triangles,
3641 const std::vector<Vec4I>& quads,
3644 return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3645 bandWidth, bandWidth,
true);
3653 inline std::ostream&
3656 ostr <<
"{[ " << rhs.
mXPrim <<
", " << rhs.
mXDist <<
"]";
3657 ostr <<
" [ " << rhs.
mYPrim <<
", " << rhs.
mYDist <<
"]";
3658 ostr <<
" [ " << rhs.
mZPrim <<
", " << rhs.
mZDist <<
"]}";
3678 const std::vector<Vec3s>& pointList,
3679 const std::vector<Vec4I>& polygonList);
3681 void run(
bool threaded =
true);
3684 inline void operator() (
const tbb::blocked_range<size_t> &range);
3692 struct Primitive {
Vec3d a, b, c, d;
Int32 index; };
3694 template<
bool IsQuad>
3695 inline void voxelize(
const Primitive&);
3697 template<
bool IsQuad>
3698 inline bool evalPrimitive(
const Coord&,
const Primitive&);
3700 inline bool rayTriangleIntersection(
const Vec3d& origin,
const Vec3d& dir,
3707 const std::vector<Vec3s>& mPointList;
3708 const std::vector<Vec4I>& mPolygonList;
3712 IntTreeT mLastPrimTree;
3718 MeshToVoxelEdgeData::GenEdgeData::GenEdgeData(
3719 const std::vector<Vec3s>& pointList,
3720 const std::vector<Vec4I>& polygonList)
3723 , mPointList(pointList)
3724 , mPolygonList(polygonList)
3726 , mLastPrimAccessor(mLastPrimTree)
3735 , mPointList(rhs.mPointList)
3736 , mPolygonList(rhs.mPolygonList)
3738 , mLastPrimAccessor(mLastPrimTree)
3747 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPolygonList.size()), *
this);
3749 (*this)(tbb::blocked_range<size_t>(0, mPolygonList.size()));
3758 using NodeChainType = RootNodeType::NodeChainType;
3759 static_assert(boost::mpl::size<NodeChainType>::value > 1,
"expected tree height > 1");
3760 using InternalNodeType = boost::mpl::at<NodeChainType, boost::mpl::int_<1> >::type;
3768 for ( ; leafIt; ++leafIt) {
3769 ijk = leafIt->origin();
3776 InternalNodeType* node = rhs.mAccessor.
getNode<InternalNodeType>();
3778 rhs.mAccessor.
clear();
3788 if (!lhsLeafPt->isValueOn(offset)) {
3789 lhsLeafPt->setValueOn(offset, rhsValue);
3821 for (
size_t n = range.begin(); n < range.end(); ++n) {
3823 const Vec4I& verts = mPolygonList[n];
3825 prim.index =
Int32(n);
3826 prim.a =
Vec3d(mPointList[verts[0]]);
3827 prim.b =
Vec3d(mPointList[verts[1]]);
3828 prim.c =
Vec3d(mPointList[verts[2]]);
3831 prim.d =
Vec3d(mPointList[verts[3]]);
3832 voxelize<true>(prim);
3834 voxelize<false>(prim);
3840 template<
bool IsQuad>
3842 MeshToVoxelEdgeData::GenEdgeData::voxelize(
const Primitive& prim)
3844 std::deque<Coord> coordList;
3847 ijk = Coord::floor(prim.a);
3848 coordList.push_back(ijk);
3850 evalPrimitive<IsQuad>(ijk, prim);
3852 while (!coordList.empty()) {
3854 ijk = coordList.back();
3855 coordList.pop_back();
3857 for (
Int32 i = 0; i < 26; ++i) {
3860 if (prim.index != mLastPrimAccessor.
getValue(nijk)) {
3861 mLastPrimAccessor.
setValue(nijk, prim.index);
3862 if(evalPrimitive<IsQuad>(nijk, prim)) coordList.push_back(nijk);
3869 template<
bool IsQuad>
3871 MeshToVoxelEdgeData::GenEdgeData::evalPrimitive(
const Coord& ijk,
const Primitive& prim)
3873 Vec3d uvw, org(ijk[0], ijk[1], ijk[2]);
3874 bool intersecting =
false;
3881 double dist = (org -
3884 if (rayTriangleIntersection(org,
Vec3d(1.0, 0.0, 0.0), prim.a, prim.c, prim.b, t)) {
3885 if (t < edgeData.
mXDist) {
3886 edgeData.
mXDist = float(t);
3887 edgeData.
mXPrim = prim.index;
3888 intersecting =
true;
3892 if (rayTriangleIntersection(org,
Vec3d(0.0, 1.0, 0.0), prim.a, prim.c, prim.b, t)) {
3893 if (t < edgeData.
mYDist) {
3894 edgeData.
mYDist = float(t);
3895 edgeData.
mYPrim = prim.index;
3896 intersecting =
true;
3900 if (rayTriangleIntersection(org,
Vec3d(0.0, 0.0, 1.0), prim.a, prim.c, prim.b, t)) {
3901 if (t < edgeData.
mZDist) {
3902 edgeData.
mZDist = float(t);
3903 edgeData.
mZPrim = prim.index;
3904 intersecting =
true;
3910 double secondDist = (org -
3913 if (secondDist < dist) dist = secondDist;
3915 if (rayTriangleIntersection(org,
Vec3d(1.0, 0.0, 0.0), prim.a, prim.d, prim.c, t)) {
3916 if (t < edgeData.
mXDist) {
3917 edgeData.
mXDist = float(t);
3918 edgeData.
mXPrim = prim.index;
3919 intersecting =
true;
3923 if (rayTriangleIntersection(org,
Vec3d(0.0, 1.0, 0.0), prim.a, prim.d, prim.c, t)) {
3924 if (t < edgeData.
mYDist) {
3925 edgeData.
mYDist = float(t);
3926 edgeData.
mYPrim = prim.index;
3927 intersecting =
true;
3931 if (rayTriangleIntersection(org,
Vec3d(0.0, 0.0, 1.0), prim.a, prim.d, prim.c, t)) {
3932 if (t < edgeData.
mZDist) {
3933 edgeData.
mZDist = float(t);
3934 edgeData.
mZPrim = prim.index;
3935 intersecting =
true;
3940 if (intersecting) mAccessor.
setValue(ijk, edgeData);
3942 return (dist < 0.86602540378443861);
3947 MeshToVoxelEdgeData::GenEdgeData::rayTriangleIntersection(
3958 double divisor = s1.
dot(e1);
3959 if (!(std::abs(divisor) > 0.0))
return false;
3963 double inv_divisor = 1.0 / divisor;
3964 Vec3d d = origin - a;
3965 double b1 = d.
dot(s1) * inv_divisor;
3967 if (b1 < 0.0 || b1 > 1.0)
return false;
3970 double b2 = dir.
dot(s2) * inv_divisor;
3972 if (b2 < 0.0 || (b1 + b2) > 1.0)
return false;
3976 t = e2.dot(s2) * inv_divisor;
3977 return (t < 0.0) ?
false :
true;
3993 const std::vector<Vec3s>& pointList,
3994 const std::vector<Vec4I>& polygonList)
4008 std::vector<Vec3d>& points,
4009 std::vector<Index32>& primitives)
4019 point[0] = double(coord[0]) + data.
mXDist;
4020 point[1] = double(coord[1]);
4021 point[2] = double(coord[2]);
4023 points.push_back(point);
4024 primitives.push_back(data.
mXPrim);
4028 point[0] = double(coord[0]);
4029 point[1] = double(coord[1]) + data.
mYDist;
4030 point[2] = double(coord[2]);
4032 points.push_back(point);
4033 primitives.push_back(data.
mYPrim);
4037 point[0] = double(coord[0]);
4038 point[1] = double(coord[1]);
4039 point[2] = double(coord[2]) + data.
mZDist;
4041 points.push_back(point);
4042 primitives.push_back(data.
mZPrim);
4052 point[0] = double(coord[0]);
4053 point[1] = double(coord[1]) + data.
mYDist;
4054 point[2] = double(coord[2]);
4056 points.push_back(point);
4057 primitives.push_back(data.
mYPrim);
4061 point[0] = double(coord[0]);
4062 point[1] = double(coord[1]);
4063 point[2] = double(coord[2]) + data.
mZDist;
4065 points.push_back(point);
4066 primitives.push_back(data.
mZPrim);
4074 point[0] = double(coord[0]);
4075 point[1] = double(coord[1]) + data.
mYDist;
4076 point[2] = double(coord[2]);
4078 points.push_back(point);
4079 primitives.push_back(data.
mYPrim);
4088 point[0] = double(coord[0]) + data.
mXDist;
4089 point[1] = double(coord[1]);
4090 point[2] = double(coord[2]);
4092 points.push_back(point);
4093 primitives.push_back(data.
mXPrim);
4097 point[0] = double(coord[0]);
4098 point[1] = double(coord[1]) + data.
mYDist;
4099 point[2] = double(coord[2]);
4101 points.push_back(point);
4102 primitives.push_back(data.
mYPrim);
4112 point[0] = double(coord[0]) + data.
mXDist;
4113 point[1] = double(coord[1]);
4114 point[2] = double(coord[2]);
4116 points.push_back(point);
4117 primitives.push_back(data.
mXPrim);
4126 point[0] = double(coord[0]) + data.
mXDist;
4127 point[1] = double(coord[1]);
4128 point[2] = double(coord[2]);
4130 points.push_back(point);
4131 primitives.push_back(data.
mXPrim);
4135 point[0] = double(coord[0]);
4136 point[1] = double(coord[1]);
4137 point[2] = double(coord[2]) + data.
mZDist;
4139 points.push_back(point);
4140 primitives.push_back(data.
mZPrim);
4149 point[0] = double(coord[0]);
4150 point[1] = double(coord[1]);
4151 point[2] = double(coord[2]) + data.
mZDist;
4153 points.push_back(point);
4154 primitives.push_back(data.
mZPrim);
4160 template<
typename Gr
idType,
typename VecType>
4161 inline typename GridType::Ptr
4163 const openvdb::math::Transform& xform,
4164 typename VecType::ValueType halfWidth)
4170 points[0] =
Vec3s(pmin[0], pmin[1], pmin[2]);
4171 points[1] =
Vec3s(pmin[0], pmin[1], pmax[2]);
4172 points[2] =
Vec3s(pmax[0], pmin[1], pmax[2]);
4173 points[3] =
Vec3s(pmax[0], pmin[1], pmin[2]);
4174 points[4] =
Vec3s(pmin[0], pmax[1], pmin[2]);
4175 points[5] =
Vec3s(pmin[0], pmax[1], pmax[2]);
4176 points[6] =
Vec3s(pmax[0], pmax[1], pmax[2]);
4177 points[7] =
Vec3s(pmax[0], pmax[1], pmin[2]);
4180 faces[0] =
Vec4I(0, 1, 2, 3);
4181 faces[1] =
Vec4I(7, 6, 5, 4);
4182 faces[2] =
Vec4I(4, 5, 1, 0);
4183 faces[3] =
Vec4I(6, 7, 3, 2);
4184 faces[4] =
Vec4I(0, 3, 7, 4);
4185 faces[5] =
Vec4I(1, 5, 6, 2);
4189 return meshToVolume<GridType>(mesh, xform, halfWidth, halfWidth);
4197 #endif // OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
Vec3< float > Vec3s
Definition: Vec3.h:707
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
uint32_t Index32
Definition: Types.h:55
LeafNodeT * probeLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains voxel (x, y, z), or nullptr 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
Tree< typename RootNodeType::template ValueConverter< Int32 >::Type > Type
Definition: Tree.h:226
#define OPENVDB_LOG_DEBUG(message)
In debug builds only, log a debugging message of the form 'someVar << "text" << ...'.
Definition: logging.h:288
Convert polygonal meshes that consist of quads and/or triangles into signed or unsigned distance fiel...
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:101
Efficient multi-threaded replacement of the background values in tree.
Axis-aligned bounding box.
Definition: BBox.h:47
OPENVDB_API const Coord COORD_OFFSETS[26]
coordinate offset table for neighboring voxels
ValueConverter<T>::Type is the type of a tree having the same hierarchy as this tree but a different ...
Definition: Tree.h:225
void setValueOnly(const Coord &xyz, const ValueType &value)
Set the value of the voxel at the given coordinate but don't change its active state.
Definition: ValueAccessor.h:296
const ValueT & getValue() const
Return the tile or voxel value to which this iterator is currently pointing.
Definition: TreeIterator.h:734
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
Defined various multi-threaded utility functions for trees.
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: ValueAccessor.h:256
static const Real LEVEL_SET_HALF_WIDTH
Definition: Types.h:270
Real GodunovsNormSqrd(bool isOutside, Real dP_xm, Real dP_xp, Real dP_ym, Real dP_yp, Real dP_zm, Real dP_zp)
Definition: FiniteDifference.h:353
Vec2< T > maxComponent(const Vec2< T > &v1, const Vec2< T > &v2)
Return component-wise maximum of the two vectors.
Definition: Vec2.h:562
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:727
#define OPENVDB_VERSION_NAME
Definition: version.h:43
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 isValueOn(const Coord &xyz) const
Return the active state of the voxel at the given coordinates.
Definition: ValueAccessor.h:263
Vec3< double > Vec3d
Definition: Vec3.h:708
Definition: Exceptions.h:91
Propagates the sign of distance values from the active voxels in the narrow band to the inactive valu...
Base class for tree-traversal iterators over all leaf nodes (but not leaf voxels) ...
Definition: TreeIterator.h:1228
const LeafNodeT * probeConstLeaf(const Coord &xyz) const
Return a pointer to the leaf node that contains voxel (x, y, z), or nullptr if no such node exists...
Definition: ValueAccessor.h:429
T dot(const Vec3< T > &v) const
Dot product.
Definition: Vec3.h:215
bool normalize(T eps=T(1.0e-7))
this = normalized this
Definition: Vec3.h:376
Definition: Exceptions.h:39
void clear()
Remove all tiles from this tree and all nodes other than the root node.
Definition: Tree.h:1489
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:76
void expand(ElementType padding)
Pad this bounding box.
Definition: BBox.h:356
typename RootNodeType::LeafNodeType LeafNodeType
Definition: Tree.h:214
void clearAllAccessors()
Clear all registered accessors.
Definition: Tree.h:1548
virtual void clear()
Remove all nodes from this cache, then reinsert the root node.
Definition: ValueAccessor.h:438
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
math::Vec4< Index32 > Vec4I
Definition: Types.h:91
Base class for tree-traversal iterators over tile and voxel values.
Definition: TreeIterator.h:658
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:407
std::shared_ptr< T > SharedPtr
Definition: Types.h:130
LeafIter beginLeaf()
Return an iterator over all leaf nodes in this tree.
Definition: Tree.h:1180
bool operator<(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition: Tuple.h:174
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:52
const Vec3T & min() const
Return a const reference to the minimum point of the BBox.
Definition: BBox.h:84
void merge(Tree &other, MergePolicy=MERGE_ACTIVE_STATES)
Efficiently merge another tree into this tree using one of several schemes.
Definition: Tree.h:1865
LeafNodeType * probeLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains voxel (x, y, z). If no such node exists, return nullptr.
Definition: Tree.h:1735
const Vec3T & max() const
Return a const reference to the maximum point of the BBox.
Definition: BBox.h:87
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
int32_t Int32
Definition: Types.h:59
Index32 Index
Definition: Types.h:57
Vec2< T > minComponent(const Vec2< T > &v1, const Vec2< T > &v2)
Return component-wise minimum of the two vectors.
Definition: Vec2.h:553
bool probeValue(const Coord &xyz, ValueType &value) const
Return the active state of the voxel as well as its value.
Definition: ValueAccessor.h:266
OPENVDB_API const Index32 INVALID_IDX
_RootNodeType RootNodeType
Definition: Tree.h:211
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:324
bool operator>(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition: Tuple.h:186
void setValue(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:287
Type Pow2(Type x)
Return .
Definition: Math.h:514
Vec3< T > cross(const Vec3< T > &v) const
Return the cross product of "this" vector and v;.
Definition: Vec3.h:244