16 #ifndef OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED 17 #define OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED 25 #include <openvdb/thread/Threading.h> 32 #include <tbb/blocked_range.h> 33 #include <tbb/enumerable_thread_specific.h> 34 #include <tbb/parallel_for.h> 35 #include <tbb/parallel_reduce.h> 36 #include <tbb/partitioner.h> 37 #include <tbb/task_group.h> 38 #include <tbb/task_arena.h> 46 #include <type_traits> 127 template <
typename Gr
idType,
typename MeshDataAdapter,
typename InteriorTest = std::
nullptr_t>
128 typename GridType::Ptr
132 float exteriorBandWidth = 3.0f,
133 float interiorBandWidth = 3.0f,
136 InteriorTest interiorTest =
nullptr,
158 template <
typename Gr
idType,
typename MeshDataAdapter,
typename Interrupter,
typename InteriorTest = std::
nullptr_t>
159 typename GridType::Ptr
161 Interrupter& interrupter,
164 float exteriorBandWidth = 3.0f,
165 float interiorBandWidth = 3.0f,
168 InteriorTest interiorTest =
nullptr,
185 template<
typename Po
intType,
typename PolygonType>
189 const std::vector<PolygonType>& polygons)
190 : mPointArray(points.empty() ? nullptr : &points[0])
191 , mPointArraySize(points.size())
192 , mPolygonArray(polygons.empty() ? nullptr : &polygons[0])
193 , mPolygonArraySize(polygons.size())
198 const PolygonType* polygonArray,
size_t polygonArraySize)
199 : mPointArray(pointArray)
200 , mPointArraySize(pointArraySize)
201 , mPolygonArray(polygonArray)
202 , mPolygonArraySize(polygonArraySize)
211 return (PolygonType::size == 3 || mPolygonArray[n][3] ==
util::INVALID_IDX) ? 3 : 4;
217 const PointType& p = mPointArray[mPolygonArray[n][int(v)]];
218 pos[0] = double(p[0]);
219 pos[1] = double(p[1]);
220 pos[2] = double(p[2]);
224 PointType
const *
const mPointArray;
225 size_t const mPointArraySize;
226 PolygonType
const *
const mPolygonArray;
227 size_t const mPolygonArraySize;
255 template<
typename Gr
idType>
256 typename GridType::Ptr
258 const openvdb::math::Transform& xform,
259 const std::vector<Vec3s>& points,
260 const std::vector<Vec3I>& triangles,
264 template<
typename Gr
idType,
typename Interrupter>
265 typename GridType::Ptr
267 Interrupter& interrupter,
268 const openvdb::math::Transform& xform,
269 const std::vector<Vec3s>& points,
270 const std::vector<Vec3I>& triangles,
289 template<
typename Gr
idType>
290 typename GridType::Ptr
292 const openvdb::math::Transform& xform,
293 const std::vector<Vec3s>& points,
294 const std::vector<Vec4I>& quads,
298 template<
typename Gr
idType,
typename Interrupter>
299 typename GridType::Ptr
301 Interrupter& interrupter,
302 const openvdb::math::Transform& xform,
303 const std::vector<Vec3s>& points,
304 const std::vector<Vec4I>& quads,
324 template<
typename Gr
idType>
325 typename GridType::Ptr
327 const openvdb::math::Transform& xform,
328 const std::vector<Vec3s>& points,
329 const std::vector<Vec3I>& triangles,
330 const std::vector<Vec4I>& quads,
334 template<
typename Gr
idType,
typename Interrupter>
335 typename GridType::Ptr
337 Interrupter& interrupter,
338 const openvdb::math::Transform& xform,
339 const std::vector<Vec3s>& points,
340 const std::vector<Vec3I>& triangles,
341 const std::vector<Vec4I>& quads,
363 template<
typename Gr
idType>
364 typename GridType::Ptr
366 const openvdb::math::Transform& xform,
367 const std::vector<Vec3s>& points,
368 const std::vector<Vec3I>& triangles,
369 const std::vector<Vec4I>& quads,
374 template<
typename Gr
idType,
typename Interrupter>
375 typename GridType::Ptr
377 Interrupter& interrupter,
378 const openvdb::math::Transform& xform,
379 const std::vector<Vec3s>& points,
380 const std::vector<Vec3I>& triangles,
381 const std::vector<Vec4I>& quads,
400 template<
typename Gr
idType>
401 typename GridType::Ptr
403 const openvdb::math::Transform& xform,
404 const std::vector<Vec3s>& points,
405 const std::vector<Vec3I>& triangles,
406 const std::vector<Vec4I>& quads,
410 template<
typename Gr
idType,
typename Interrupter>
411 typename GridType::Ptr
413 Interrupter& interrupter,
414 const openvdb::math::Transform& xform,
415 const std::vector<Vec3s>& points,
416 const std::vector<Vec3I>& triangles,
417 const std::vector<Vec4I>& quads,
430 template<
typename Gr
idType,
typename VecType>
431 typename GridType::Ptr
433 const openvdb::math::Transform& xform,
446 template <
typename FloatTreeT>
464 : mXDist(dist), mYDist(dist), mZDist(dist)
507 void convert(
const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList);
512 void getEdgeData(
Accessor& acc,
const Coord& ijk,
513 std::vector<Vec3d>& points, std::vector<Index32>& primitives);
534 namespace mesh_to_volume_internal {
536 template<
typename Po
intType>
537 struct TransformPoints {
539 TransformPoints(
const PointType* pointsIn, PointType* pointsOut,
541 : mPointsIn(pointsIn), mPointsOut(pointsOut), mXform(&xform)
545 void operator()(
const tbb::blocked_range<size_t>& range)
const {
549 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
551 const PointType& wsP = mPointsIn[n];
552 pos[0] = double(wsP[0]);
553 pos[1] = double(wsP[1]);
554 pos[2] = double(wsP[2]);
556 pos = mXform->worldToIndex(pos);
558 PointType& isP = mPointsOut[n];
559 isP[0] =
typename PointType::value_type(pos[0]);
560 isP[1] =
typename PointType::value_type(pos[1]);
561 isP[2] =
typename PointType::value_type(pos[2]);
565 PointType
const *
const mPointsIn;
566 PointType *
const mPointsOut;
571 template<
typename ValueType>
574 static ValueType epsilon() {
return ValueType(1e-7); }
575 static ValueType minNarrowBandWidth() {
return ValueType(1.0 + 1e-6); }
582 template<
typename TreeType>
583 class CombineLeafNodes
590 using Int32LeafNodeType =
typename Int32TreeType::LeafNodeType;
592 CombineLeafNodes(
TreeType& lhsDistTree, Int32TreeType& lhsIdxTree,
593 LeafNodeType ** rhsDistNodes, Int32LeafNodeType ** rhsIdxNodes)
594 : mDistTree(&lhsDistTree)
595 , mIdxTree(&lhsIdxTree)
596 , mRhsDistNodes(rhsDistNodes)
597 , mRhsIdxNodes(rhsIdxNodes)
601 void operator()(
const tbb::blocked_range<size_t>& range)
const {
606 using DistValueType =
typename LeafNodeType::ValueType;
607 using IndexValueType =
typename Int32LeafNodeType::ValueType;
609 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
611 const Coord& origin = mRhsDistNodes[n]->origin();
613 LeafNodeType* lhsDistNode = distAcc.
probeLeaf(origin);
614 Int32LeafNodeType* lhsIdxNode = idxAcc.
probeLeaf(origin);
616 DistValueType* lhsDistData = lhsDistNode->buffer().data();
617 IndexValueType* lhsIdxData = lhsIdxNode->buffer().data();
619 const DistValueType* rhsDistData = mRhsDistNodes[n]->buffer().data();
620 const IndexValueType* rhsIdxData = mRhsIdxNodes[n]->buffer().data();
623 for (
Index32 offset = 0; offset < LeafNodeType::SIZE; ++offset) {
627 const DistValueType& lhsValue = lhsDistData[offset];
628 const DistValueType& rhsValue = rhsDistData[offset];
630 if (rhsValue < lhsValue) {
631 lhsDistNode->setValueOn(offset, rhsValue);
632 lhsIdxNode->setValueOn(offset, rhsIdxData[offset]);
634 lhsIdxNode->setValueOn(offset,
635 std::min(lhsIdxData[offset], rhsIdxData[offset]));
640 delete mRhsDistNodes[n];
641 delete mRhsIdxNodes[n];
648 Int32TreeType *
const mIdxTree;
650 LeafNodeType **
const mRhsDistNodes;
651 Int32LeafNodeType **
const mRhsIdxNodes;
658 template<
typename TreeType>
659 struct StashOriginAndStoreOffset
663 StashOriginAndStoreOffset(std::vector<LeafNodeType*>& nodes, Coord* coordinates)
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 mCoordinates[n] = origin;
672 origin[0] =
static_cast<int>(n);
676 LeafNodeType **
const mNodes;
677 Coord *
const mCoordinates;
681 template<
typename TreeType>
686 RestoreOrigin(std::vector<LeafNodeType*>& nodes,
const Coord* coordinates)
687 : mNodes(nodes.empty() ?
nullptr : &nodes[0]), mCoordinates(coordinates)
691 void operator()(
const tbb::blocked_range<size_t>& range)
const {
692 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
693 Coord& origin =
const_cast<Coord&
>(mNodes[n]->origin());
694 origin[0] = mCoordinates[n][0];
698 LeafNodeType **
const mNodes;
699 Coord
const *
const mCoordinates;
703 template<
typename TreeType>
704 class ComputeNodeConnectivity
709 ComputeNodeConnectivity(
const TreeType& tree,
const Coord* coordinates,
710 size_t* offsets,
size_t numNodes,
const CoordBBox& bbox)
712 , mCoordinates(coordinates)
714 , mNumNodes(numNodes)
719 ComputeNodeConnectivity(
const ComputeNodeConnectivity&) =
default;
722 ComputeNodeConnectivity& operator=(
const ComputeNodeConnectivity&) =
delete;
724 void operator()(
const tbb::blocked_range<size_t>& range)
const {
726 size_t* offsetsNextX = mOffsets;
727 size_t* offsetsPrevX = mOffsets + mNumNodes;
728 size_t* offsetsNextY = mOffsets + mNumNodes * 2;
729 size_t* offsetsPrevY = mOffsets + mNumNodes * 3;
730 size_t* offsetsNextZ = mOffsets + mNumNodes * 4;
731 size_t* offsetsPrevZ = mOffsets + mNumNodes * 5;
735 const Int32 DIM =
static_cast<Int32>(LeafNodeType::DIM);
737 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
738 const Coord& origin = mCoordinates[n];
739 offsetsNextX[n] = findNeighbourNode(acc, origin, Coord(DIM, 0, 0));
740 offsetsPrevX[n] = findNeighbourNode(acc, origin, Coord(-DIM, 0, 0));
741 offsetsNextY[n] = findNeighbourNode(acc, origin, Coord(0, DIM, 0));
742 offsetsPrevY[n] = findNeighbourNode(acc, origin, Coord(0, -DIM, 0));
743 offsetsNextZ[n] = findNeighbourNode(acc, origin, Coord(0, 0, DIM));
744 offsetsPrevZ[n] = findNeighbourNode(acc, origin, Coord(0, 0, -DIM));
749 const Coord& start,
const Coord& step)
const 751 Coord ijk = start + step;
754 while (bbox.isInside(ijk)) {
756 if (node)
return static_cast<size_t>(node->origin()[0]);
766 Coord
const *
const mCoordinates;
767 size_t *
const mOffsets;
769 const size_t mNumNodes;
774 template<
typename TreeType>
775 struct LeafNodeConnectivityTable
781 LeafNodeConnectivityTable(
TreeType& tree)
786 if (mLeafNodes.empty())
return;
791 const tbb::blocked_range<size_t> range(0, mLeafNodes.size());
795 std::unique_ptr<Coord[]> coordinates{
new Coord[mLeafNodes.size()]};
796 tbb::parallel_for(range,
797 StashOriginAndStoreOffset<TreeType>(mLeafNodes, coordinates.get()));
800 mOffsets.reset(
new size_t[mLeafNodes.size() * 6]);
803 tbb::parallel_for(range, ComputeNodeConnectivity<TreeType>(
804 tree, coordinates.get(), mOffsets.get(), mLeafNodes.size(), bbox));
807 tbb::parallel_for(range, RestoreOrigin<TreeType>(mLeafNodes, coordinates.get()));
810 size_t size()
const {
return mLeafNodes.size(); }
812 std::vector<LeafNodeType*>& nodes() {
return mLeafNodes; }
813 const std::vector<LeafNodeType*>& nodes()
const {
return mLeafNodes; }
816 const size_t* offsetsNextX()
const {
return mOffsets.get(); }
817 const size_t* offsetsPrevX()
const {
return mOffsets.get() + mLeafNodes.size(); }
819 const size_t* offsetsNextY()
const {
return mOffsets.get() + mLeafNodes.size() * 2; }
820 const size_t* offsetsPrevY()
const {
return mOffsets.get() + mLeafNodes.size() * 3; }
822 const size_t* offsetsNextZ()
const {
return mOffsets.get() + mLeafNodes.size() * 4; }
823 const size_t* offsetsPrevZ()
const {
return mOffsets.get() + mLeafNodes.size() * 5; }
826 std::vector<LeafNodeType*> mLeafNodes;
827 std::unique_ptr<size_t[]> mOffsets;
831 template<
typename TreeType>
832 class SweepExteriorSign
840 using ConnectivityTable = LeafNodeConnectivityTable<TreeType>;
842 SweepExteriorSign(
Axis axis,
const std::vector<size_t>& startNodeIndices,
843 ConnectivityTable& connectivity)
844 : mStartNodeIndices(startNodeIndices.empty() ?
nullptr : &startNodeIndices[0])
845 , mConnectivity(&connectivity)
850 void operator()(
const tbb::blocked_range<size_t>& range)
const {
852 constexpr
Int32 DIM =
static_cast<Int32>(LeafNodeType::DIM);
854 std::vector<LeafNodeType*>& nodes = mConnectivity->nodes();
857 size_t idxA = 0, idxB = 1;
860 const size_t* nextOffsets = mConnectivity->offsetsNextZ();
861 const size_t* prevOffsets = mConnectivity->offsetsPrevZ();
869 nextOffsets = mConnectivity->offsetsNextY();
870 prevOffsets = mConnectivity->offsetsPrevY();
872 }
else if (mAxis ==
X_AXIS) {
878 nextOffsets = mConnectivity->offsetsNextX();
879 prevOffsets = mConnectivity->offsetsPrevX();
884 Int32& a = ijk[idxA];
885 Int32& b = ijk[idxB];
887 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
889 size_t startOffset = mStartNodeIndices[n];
890 size_t lastOffset = startOffset;
894 for (a = 0; a < DIM; ++a) {
895 for (b = 0; b < DIM; ++b) {
897 pos =
static_cast<Int32>(LeafNodeType::coordToOffset(ijk));
898 size_t offset = startOffset;
901 while ( offset != ConnectivityTable::INVALID_OFFSET &&
902 traceVoxelLine(*nodes[offset], pos, step) ) {
905 offset = nextOffsets[offset];
910 while (offset != ConnectivityTable::INVALID_OFFSET) {
912 offset = nextOffsets[offset];
917 pos += step * (DIM - 1);
918 while ( offset != ConnectivityTable::INVALID_OFFSET &&
919 traceVoxelLine(*nodes[offset], pos, -step)) {
920 offset = prevOffsets[offset];
928 bool traceVoxelLine(LeafNodeType& node,
Int32 pos,
const Int32 step)
const {
930 ValueType* data = node.buffer().data();
932 bool isOutside =
true;
934 for (
Index i = 0; i < LeafNodeType::DIM; ++i) {
937 ValueType& dist = data[pos];
939 if (dist < ValueType(0.0)) {
943 if (!(dist > ValueType(0.75))) isOutside =
false;
945 if (isOutside) dist = ValueType(-dist);
956 size_t const *
const mStartNodeIndices;
957 ConnectivityTable *
const mConnectivity;
963 template<
typename LeafNodeType>
965 seedFill(LeafNodeType& node)
967 using ValueType =
typename LeafNodeType::ValueType;
968 using Queue = std::deque<Index>;
971 ValueType* data = node.buffer().data();
975 for (
Index pos = 0; pos < LeafNodeType::SIZE; ++pos) {
976 if (data[pos] < 0.0) seedPoints.push_back(pos);
979 if (seedPoints.empty())
return;
982 for (Queue::iterator it = seedPoints.begin(); it != seedPoints.end(); ++it) {
983 ValueType& dist = data[*it];
990 Index pos(0), nextPos(0);
992 while (!seedPoints.empty()) {
994 pos = seedPoints.back();
995 seedPoints.pop_back();
997 ValueType& dist = data[pos];
999 if (!(dist < ValueType(0.0))) {
1003 ijk = LeafNodeType::offsetToLocalCoord(pos);
1006 nextPos = pos - LeafNodeType::DIM * LeafNodeType::DIM;
1007 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1010 if (ijk[0] != (LeafNodeType::DIM - 1)) {
1011 nextPos = pos + LeafNodeType::DIM * LeafNodeType::DIM;
1012 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1016 nextPos = pos - LeafNodeType::DIM;
1017 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1020 if (ijk[1] != (LeafNodeType::DIM - 1)) {
1021 nextPos = pos + LeafNodeType::DIM;
1022 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1027 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1030 if (ijk[2] != (LeafNodeType::DIM - 1)) {
1032 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1039 template<
typename LeafNodeType>
1041 scanFill(LeafNodeType& node)
1043 bool updatedNode =
false;
1045 using ValueType =
typename LeafNodeType::ValueType;
1046 ValueType* data = node.buffer().data();
1050 bool updatedSign =
true;
1051 while (updatedSign) {
1053 updatedSign =
false;
1055 for (
Index pos = 0; pos < LeafNodeType::SIZE; ++pos) {
1057 ValueType& dist = data[pos];
1059 if (!(dist < ValueType(0.0)) && dist > ValueType(0.75)) {
1061 ijk = LeafNodeType::offsetToLocalCoord(pos);
1064 if (ijk[2] != 0 && data[pos - 1] < ValueType(0.0)) {
1066 dist = ValueType(-dist);
1069 }
else if (ijk[2] != (LeafNodeType::DIM - 1) && data[pos + 1] < ValueType(0.0)) {
1071 dist = ValueType(-dist);
1074 }
else if (ijk[1] != 0 && data[pos - LeafNodeType::DIM] < ValueType(0.0)) {
1076 dist = ValueType(-dist);
1079 }
else if (ijk[1] != (LeafNodeType::DIM - 1)
1080 && data[pos + LeafNodeType::DIM] < ValueType(0.0))
1083 dist = ValueType(-dist);
1086 }
else if (ijk[0] != 0
1087 && data[pos - LeafNodeType::DIM * LeafNodeType::DIM] < ValueType(0.0))
1090 dist = ValueType(-dist);
1093 }
else if (ijk[0] != (LeafNodeType::DIM - 1)
1094 && data[pos + LeafNodeType::DIM * LeafNodeType::DIM] < ValueType(0.0))
1097 dist = ValueType(-dist);
1102 updatedNode |= updatedSign;
1109 template<
typename TreeType>
1110 class SeedFillExteriorSign
1116 SeedFillExteriorSign(std::vector<LeafNodeType*>& nodes,
const bool* changedNodeMask)
1117 : mNodes(nodes.empty() ?
nullptr : &nodes[0])
1118 , mChangedNodeMask(changedNodeMask)
1122 void operator()(
const tbb::blocked_range<size_t>& range)
const {
1123 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1124 if (mChangedNodeMask[n]) {
1130 scanFill(*mNodes[n]);
1135 LeafNodeType **
const mNodes;
1136 const bool *
const mChangedNodeMask;
1140 template<
typename ValueType>
1143 FillArray(ValueType* array,
const ValueType v) : mArray(array), mValue(v) { }
1145 void operator()(
const tbb::blocked_range<size_t>& range)
const {
1146 const ValueType v = mValue;
1147 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1152 ValueType *
const mArray;
1153 const ValueType mValue;
1157 template<
typename ValueType>
1159 fillArray(ValueType* array,
const ValueType val,
const size_t length)
1161 const auto grainSize = std::max<size_t>(
1162 length / tbb::this_task_arena::max_concurrency(), 1024);
1163 const tbb::blocked_range<size_t> range(0, length, grainSize);
1164 tbb::parallel_for(range, FillArray<ValueType>(array, val), tbb::simple_partitioner());
1168 template<
typename TreeType>
1175 SyncVoxelMask(std::vector<LeafNodeType*>& nodes,
1176 const bool* changedNodeMask,
bool* changedVoxelMask)
1177 : mNodes(nodes.empty() ?
nullptr : &nodes[0])
1178 , mChangedNodeMask(changedNodeMask)
1179 , mChangedVoxelMask(changedVoxelMask)
1183 void operator()(
const tbb::blocked_range<size_t>& range)
const {
1184 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1186 if (mChangedNodeMask[n]) {
1187 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1189 ValueType* data = mNodes[n]->buffer().data();
1191 for (
Index pos = 0; pos < LeafNodeType::SIZE; ++pos) {
1193 data[pos] = ValueType(-data[pos]);
1201 LeafNodeType **
const mNodes;
1202 bool const *
const mChangedNodeMask;
1203 bool *
const mChangedVoxelMask;
1207 template<
typename TreeType>
1213 using ConnectivityTable = LeafNodeConnectivityTable<TreeType>;
1215 SeedPoints(ConnectivityTable& connectivity,
1216 bool* changedNodeMask,
bool* nodeMask,
bool* changedVoxelMask)
1217 : mConnectivity(&connectivity)
1218 , mChangedNodeMask(changedNodeMask)
1219 , mNodeMask(nodeMask)
1220 , mChangedVoxelMask(changedVoxelMask)
1224 void operator()(
const tbb::blocked_range<size_t>& range)
const {
1226 for (
size_t n = range.begin(), N = range.end(); n < N; ++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;
1243 bool processZ(
const size_t n,
bool firstFace)
const 1245 const size_t offset =
1246 firstFace ? mConnectivity->offsetsPrevZ()[n] : mConnectivity->offsetsNextZ()[n];
1247 if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1249 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1251 const ValueType* lhsData = mConnectivity->nodes()[n]->buffer().data();
1252 const ValueType* rhsData = mConnectivity->nodes()[offset]->buffer().data();
1254 const Index lastOffset = LeafNodeType::DIM - 1;
1255 const Index lhsOffset =
1256 firstFace ? 0 : lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1258 Index tmpPos(0), pos(0);
1259 bool changedValue =
false;
1261 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
1262 tmpPos = x << (2 * LeafNodeType::LOG2DIM);
1263 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
1264 pos = tmpPos + (y << LeafNodeType::LOG2DIM);
1266 if (lhsData[pos + lhsOffset] > ValueType(0.75)) {
1267 if (rhsData[pos + rhsOffset] < ValueType(0.0)) {
1268 changedValue =
true;
1269 mask[pos + lhsOffset] =
true;
1275 return changedValue;
1281 bool processY(
const size_t n,
bool firstFace)
const 1283 const size_t offset =
1284 firstFace ? mConnectivity->offsetsPrevY()[n] : mConnectivity->offsetsNextY()[n];
1285 if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1287 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1289 const ValueType* lhsData = mConnectivity->nodes()[n]->buffer().data();
1290 const ValueType* rhsData = mConnectivity->nodes()[offset]->buffer().data();
1292 const Index lastOffset = LeafNodeType::DIM * (LeafNodeType::DIM - 1);
1293 const Index lhsOffset =
1294 firstFace ? 0 : lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1296 Index tmpPos(0), pos(0);
1297 bool changedValue =
false;
1299 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
1300 tmpPos = x << (2 * LeafNodeType::LOG2DIM);
1301 for (
Index z = 0; z < LeafNodeType::DIM; ++z) {
1304 if (lhsData[pos + lhsOffset] > ValueType(0.75)) {
1305 if (rhsData[pos + rhsOffset] < ValueType(0.0)) {
1306 changedValue =
true;
1307 mask[pos + lhsOffset] =
true;
1313 return changedValue;
1319 bool processX(
const size_t n,
bool firstFace)
const 1321 const size_t offset =
1322 firstFace ? mConnectivity->offsetsPrevX()[n] : mConnectivity->offsetsNextX()[n];
1323 if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1325 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1327 const ValueType* lhsData = mConnectivity->nodes()[n]->buffer().data();
1328 const ValueType* rhsData = mConnectivity->nodes()[offset]->buffer().data();
1330 const Index lastOffset = LeafNodeType::DIM * LeafNodeType::DIM * (LeafNodeType::DIM-1);
1331 const Index lhsOffset =
1332 firstFace ? 0 : lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1334 Index tmpPos(0), pos(0);
1335 bool changedValue =
false;
1337 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
1338 tmpPos = y << LeafNodeType::LOG2DIM;
1339 for (
Index z = 0; z < LeafNodeType::DIM; ++z) {
1342 if (lhsData[pos + lhsOffset] > ValueType(0.75)) {
1343 if (rhsData[pos + rhsOffset] < ValueType(0.0)) {
1344 changedValue =
true;
1345 mask[pos + lhsOffset] =
true;
1351 return changedValue;
1357 ConnectivityTable *
const mConnectivity;
1358 bool *
const mChangedNodeMask;
1359 bool *
const mNodeMask;
1360 bool *
const mChangedVoxelMask;
1366 template<
typename TreeType,
typename MeshDataAdapter>
1367 struct ComputeIntersectingVoxelSign
1372 using Int32LeafNodeType =
typename Int32TreeType::LeafNodeType;
1375 using MaskArray = std::unique_ptr<bool[]>;
1376 using LocalData = std::pair<PointArray, MaskArray>;
1377 using LocalDataTable = tbb::enumerable_thread_specific<LocalData>;
1379 ComputeIntersectingVoxelSign(
1380 std::vector<LeafNodeType*>& distNodes,
1382 const Int32TreeType& indexTree,
1384 : mDistNodes(distNodes.empty() ?
nullptr : &distNodes[0])
1385 , mDistTree(&distTree)
1386 , mIndexTree(&indexTree)
1388 , mLocalDataTable(
new LocalDataTable())
1393 void operator()(
const tbb::blocked_range<size_t>& range)
const {
1400 Index xPos(0), yPos(0);
1401 Coord ijk, nijk, nodeMin, nodeMax;
1402 Vec3d cp, xyz, nxyz, dir1, dir2;
1404 LocalData& localData = mLocalDataTable->local();
1407 if (!points) points.reset(
new Vec3d[LeafNodeType::SIZE * 2]);
1409 MaskArray& mask = localData.second;
1410 if (!mask) mask.reset(
new bool[LeafNodeType::SIZE]);
1413 typename LeafNodeType::ValueOnCIter it;
1415 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1417 LeafNodeType& node = *mDistNodes[n];
1418 ValueType* data = node.buffer().data();
1420 const Int32LeafNodeType* idxNode = idxAcc.
probeConstLeaf(node.origin());
1421 const Int32* idxData = idxNode->buffer().data();
1423 nodeMin = node.origin();
1424 nodeMax = nodeMin.offsetBy(LeafNodeType::DIM - 1);
1427 memset(mask.get(), 0,
sizeof(bool) * LeafNodeType::SIZE);
1429 for (it = node.cbeginValueOn(); it; ++it) {
1430 Index pos = it.pos();
1432 ValueType& dist = data[pos];
1433 if (dist < 0.0 || dist > 0.75)
continue;
1435 ijk = node.offsetToGlobalCoord(pos);
1437 xyz[0] = double(ijk[0]);
1438 xyz[1] = double(ijk[1]);
1439 xyz[2] = double(ijk[2]);
1445 bool flipSign =
false;
1447 for (nijk[0] = bbox.min()[0]; nijk[0] <= bbox.max()[0] && !flipSign; ++nijk[0]) {
1448 xPos = (nijk[0] & (LeafNodeType::DIM - 1u)) << (2 * LeafNodeType::LOG2DIM);
1449 for (nijk[1]=bbox.min()[1]; nijk[1] <= bbox.max()[1] && !flipSign; ++nijk[1]) {
1450 yPos = xPos + ((nijk[1] & (LeafNodeType::DIM-1u)) << LeafNodeType::LOG2DIM);
1451 for (nijk[2] = bbox.min()[2]; nijk[2] <= bbox.max()[2]; ++nijk[2]) {
1452 pos = yPos + (nijk[2] & (LeafNodeType::DIM - 1u));
1454 const Int32& polyIdx = idxData[pos];
1459 const Index pointIndex = pos * 2;
1465 nxyz[0] = double(nijk[0]);
1466 nxyz[1] = double(nijk[1]);
1467 nxyz[2] = double(nijk[2]);
1469 Vec3d& point = points[pointIndex];
1471 point = closestPoint(nxyz, polyIdx);
1473 Vec3d& direction = points[pointIndex + 1];
1474 direction = nxyz - point;
1475 direction.normalize();
1478 dir1 = xyz - points[pointIndex];
1481 if (points[pointIndex + 1].dot(dir1) > 0.0) {
1492 for (
Int32 m = 0; m < 26; ++m) {
1495 if (!bbox.isInside(nijk) && distAcc.
probeValue(nijk, nval) && nval<-0.75) {
1496 nxyz[0] = double(nijk[0]);
1497 nxyz[1] = double(nijk[1]);
1498 nxyz[2] = double(nijk[2]);
1500 cp = closestPoint(nxyz, idxAcc.
getValue(nijk));
1508 if (dir2.dot(dir1) > 0.0) {
1522 Vec3d closestPoint(
const Vec3d& center,
Int32 polyIdx)
const 1524 Vec3d a, b, c, cp, uvw;
1526 const size_t polygon = size_t(polyIdx);
1527 mMesh->getIndexSpacePoint(polygon, 0, a);
1528 mMesh->getIndexSpacePoint(polygon, 1, b);
1529 mMesh->getIndexSpacePoint(polygon, 2, c);
1533 if (4 == mMesh->vertexCount(polygon)) {
1535 mMesh->getIndexSpacePoint(polygon, 3, b);
1539 if ((center - c).lengthSqr() < (center - cp).lengthSqr()) {
1548 LeafNodeType **
const mDistNodes;
1550 Int32TreeType
const *
const mIndexTree;
1560 template<
typename LeafNodeType>
1562 maskNodeInternalNeighbours(
const Index pos,
bool (&mask)[26])
1564 using NodeT = LeafNodeType;
1566 const Coord ijk = NodeT::offsetToLocalCoord(pos);
1570 mask[0] = ijk[0] != (NodeT::DIM - 1);
1572 mask[1] = ijk[0] != 0;
1574 mask[2] = ijk[1] != (NodeT::DIM - 1);
1576 mask[3] = ijk[1] != 0;
1578 mask[4] = ijk[2] != (NodeT::DIM - 1);
1580 mask[5] = ijk[2] != 0;
1584 mask[6] = mask[0] && mask[5];
1586 mask[7] = mask[1] && mask[5];
1588 mask[8] = mask[0] && mask[4];
1590 mask[9] = mask[1] && mask[4];
1592 mask[10] = mask[0] && mask[2];
1594 mask[11] = mask[1] && mask[2];
1596 mask[12] = mask[0] && mask[3];
1598 mask[13] = mask[1] && mask[3];
1600 mask[14] = mask[3] && mask[4];
1602 mask[15] = mask[3] && mask[5];
1604 mask[16] = mask[2] && mask[4];
1606 mask[17] = mask[2] && mask[5];
1610 mask[18] = mask[1] && mask[3] && mask[5];
1612 mask[19] = mask[1] && mask[3] && mask[4];
1614 mask[20] = mask[0] && mask[3] && mask[4];
1616 mask[21] = mask[0] && mask[3] && mask[5];
1618 mask[22] = mask[1] && mask[2] && mask[5];
1620 mask[23] = mask[1] && mask[2] && mask[4];
1622 mask[24] = mask[0] && mask[2] && mask[4];
1624 mask[25] = mask[0] && mask[2] && mask[5];
1628 template<
typename Compare,
typename LeafNodeType>
1630 checkNeighbours(
const Index pos,
const typename LeafNodeType::ValueType * data,
bool (&mask)[26])
1632 using NodeT = LeafNodeType;
1635 if (mask[5] && Compare::check(data[pos - 1]))
return true;
1637 if (mask[4] && Compare::check(data[pos + 1]))
return true;
1639 if (mask[3] && Compare::check(data[pos - NodeT::DIM]))
return true;
1641 if (mask[2] && Compare::check(data[pos + NodeT::DIM]))
return true;
1643 if (mask[1] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM]))
return true;
1645 if (mask[0] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM]))
return true;
1647 if (mask[6] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM]))
return true;
1649 if (mask[7] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - 1]))
return true;
1651 if (mask[8] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + 1]))
return true;
1653 if (mask[9] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + 1]))
return true;
1655 if (mask[10] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM]))
return true;
1657 if (mask[11] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM]))
return true;
1659 if (mask[12] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM]))
return true;
1661 if (mask[13] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM]))
return true;
1663 if (mask[14] && Compare::check(data[pos - NodeT::DIM + 1]))
return true;
1665 if (mask[15] && Compare::check(data[pos - NodeT::DIM - 1]))
return true;
1667 if (mask[16] && Compare::check(data[pos + NodeT::DIM + 1]))
return true;
1669 if (mask[17] && Compare::check(data[pos + NodeT::DIM - 1]))
return true;
1671 if (mask[18] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM - 1]))
return true;
1673 if (mask[19] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM + 1]))
return true;
1675 if (mask[20] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM + 1]))
return true;
1677 if (mask[21] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM - 1]))
return true;
1679 if (mask[22] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM - 1]))
return true;
1681 if (mask[23] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM + 1]))
return true;
1683 if (mask[24] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM + 1]))
return true;
1685 if (mask[25] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM - 1]))
return true;
1691 template<
typename Compare,
typename AccessorType>
1693 checkNeighbours(
const Coord& ijk, AccessorType& acc,
bool (&mask)[26])
1695 for (
Int32 m = 0; m < 26; ++m) {
1705 template<
typename TreeType>
1706 struct ValidateIntersectingVoxels
1711 struct IsNegative {
static bool check(
const ValueType v) {
return v < ValueType(0.0); } };
1713 ValidateIntersectingVoxels(
TreeType& tree, std::vector<LeafNodeType*>& nodes)
1715 , mNodes(nodes.empty() ?
nullptr : &nodes[0])
1719 void operator()(
const tbb::blocked_range<size_t>& range)
const 1722 bool neighbourMask[26];
1724 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1726 LeafNodeType& node = *mNodes[n];
1727 ValueType* data = node.buffer().data();
1729 typename LeafNodeType::ValueOnCIter it;
1730 for (it = node.cbeginValueOn(); it; ++it) {
1732 const Index pos = it.pos();
1734 ValueType& dist = data[pos];
1735 if (dist < 0.0 || dist > 0.75)
continue;
1738 maskNodeInternalNeighbours<LeafNodeType>(pos, neighbourMask);
1740 const bool hasNegativeNeighbour =
1741 checkNeighbours<IsNegative, LeafNodeType>(pos, data, neighbourMask) ||
1742 checkNeighbours<IsNegative>(node.offsetToGlobalCoord(pos), acc, neighbourMask);
1744 if (!hasNegativeNeighbour) {
1746 dist = ValueType(0.75) + Tolerance<ValueType>::epsilon();
1753 LeafNodeType **
const mNodes;
1757 template<
typename TreeType>
1758 struct RemoveSelfIntersectingSurface
1764 struct Comp {
static bool check(
const ValueType v) {
return !(v > ValueType(0.75)); } };
1766 RemoveSelfIntersectingSurface(std::vector<LeafNodeType*>& nodes,
1767 TreeType& distTree, Int32TreeType& indexTree)
1768 : mNodes(nodes.empty() ?
nullptr : &nodes[0])
1769 , mDistTree(&distTree)
1770 , mIndexTree(&indexTree)
1774 void operator()(
const tbb::blocked_range<size_t>& range)
const 1778 bool neighbourMask[26];
1780 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1782 LeafNodeType& distNode = *mNodes[n];
1783 ValueType* data = distNode.buffer().data();
1785 typename Int32TreeType::LeafNodeType* idxNode =
1788 typename LeafNodeType::ValueOnCIter it;
1789 for (it = distNode.cbeginValueOn(); it; ++it) {
1791 const Index pos = it.pos();
1793 if (!(data[pos] > 0.75))
continue;
1796 maskNodeInternalNeighbours<LeafNodeType>(pos, neighbourMask);
1798 const bool hasBoundaryNeighbour =
1799 checkNeighbours<Comp, LeafNodeType>(pos, data, neighbourMask) ||
1800 checkNeighbours<Comp>(distNode.offsetToGlobalCoord(pos),distAcc,neighbourMask);
1802 if (!hasBoundaryNeighbour) {
1803 distNode.setValueOff(pos);
1804 idxNode->setValueOff(pos);
1810 LeafNodeType * *
const mNodes;
1812 Int32TreeType *
const mIndexTree;
1819 template<
typename NodeType>
1820 struct ReleaseChildNodes
1822 ReleaseChildNodes(NodeType ** nodes) : mNodes(nodes) {}
1824 void operator()(
const tbb::blocked_range<size_t>& range)
const {
1826 using NodeMaskType =
typename NodeType::NodeMaskType;
1828 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1829 const_cast<NodeMaskType&
>(mNodes[n]->getChildMask()).setOff();
1833 NodeType **
const mNodes;
1837 template<
typename TreeType>
1842 using NodeChainType =
typename RootNodeType::NodeChainType;
1843 using InternalNodeType =
typename NodeChainType::template Get<1>;
1845 std::vector<InternalNodeType*> nodes;
1848 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
1849 ReleaseChildNodes<InternalNodeType>(nodes.empty() ?
nullptr : &nodes[0]));
1853 template<
typename TreeType>
1854 struct StealUniqueLeafNodes
1859 std::vector<LeafNodeType*>& overlappingNodes)
1860 : mLhsTree(&lhsTree)
1861 , mRhsTree(&rhsTree)
1862 , mNodes(&overlappingNodes)
1866 void operator()()
const {
1868 std::vector<LeafNodeType*> rhsLeafNodes;
1870 rhsLeafNodes.reserve(mRhsTree->leafCount());
1873 mRhsTree->stealNodes(rhsLeafNodes);
1877 for (
size_t n = 0, N = rhsLeafNodes.size(); n < N; ++n) {
1878 if (!acc.
probeLeaf(rhsLeafNodes[n]->origin())) {
1881 mNodes->push_back(rhsLeafNodes[n]);
1889 std::vector<LeafNodeType*> *
const mNodes;
1893 template<
typename DistTreeType,
typename IndexTreeType>
1895 combineData(DistTreeType& lhsDist, IndexTreeType& lhsIdx,
1896 DistTreeType& rhsDist, IndexTreeType& rhsIdx)
1898 using DistLeafNodeType =
typename DistTreeType::LeafNodeType;
1899 using IndexLeafNodeType =
typename IndexTreeType::LeafNodeType;
1901 std::vector<DistLeafNodeType*> overlappingDistNodes;
1902 std::vector<IndexLeafNodeType*> overlappingIdxNodes;
1905 tbb::task_group tasks;
1906 tasks.run(StealUniqueLeafNodes<DistTreeType>(lhsDist, rhsDist, overlappingDistNodes));
1907 tasks.run(StealUniqueLeafNodes<IndexTreeType>(lhsIdx, rhsIdx, overlappingIdxNodes));
1911 if (!overlappingDistNodes.empty() && !overlappingIdxNodes.empty()) {
1912 tbb::parallel_for(tbb::blocked_range<size_t>(0, overlappingDistNodes.size()),
1913 CombineLeafNodes<DistTreeType>(lhsDist, lhsIdx,
1914 &overlappingDistNodes[0], &overlappingIdxNodes[0]));
1924 template<
typename TreeType>
1925 struct VoxelizationData {
1927 using Ptr = std::unique_ptr<VoxelizationData>;
1931 using UCharTreeType =
typename TreeType::template ValueConverter<unsigned char>::Type;
1942 , indexAcc(indexTree)
1943 , primIdTree(MaxPrimId)
1944 , primIdAcc(primIdTree)
1950 FloatTreeAcc distAcc;
1952 Int32TreeType indexTree;
1953 Int32TreeAcc indexAcc;
1955 UCharTreeType primIdTree;
1956 UCharTreeAcc primIdAcc;
1958 unsigned char getNewPrimId() {
1974 if (mPrimCount == MaxPrimId || primIdTree.leafCount() > 1000) {
1976 primIdTree.
root().clear();
1977 primIdTree.clearAllAccessors();
1978 assert(mPrimCount == 0);
1981 return mPrimCount++;
1986 enum { MaxPrimId = 100 };
1988 unsigned char mPrimCount;
1992 template<
typename TreeType,
typename MeshDataAdapter,
typename Interrupter = util::NullInterrupter>
1993 class VoxelizePolygons
1997 using VoxelizationDataType = VoxelizationData<TreeType>;
1998 using DataTable = tbb::enumerable_thread_specific<typename VoxelizationDataType::Ptr>;
2000 VoxelizePolygons(DataTable& dataTable,
2002 Interrupter* interrupter =
nullptr)
2003 : mDataTable(&dataTable)
2005 , mInterrupter(interrupter)
2009 void operator()(
const tbb::blocked_range<size_t>& range)
const {
2011 typename VoxelizationDataType::Ptr& dataPtr = mDataTable->local();
2012 if (!dataPtr) dataPtr.reset(
new VoxelizationDataType());
2016 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2019 thread::cancelGroupExecution();
2023 const size_t numVerts = mMesh->vertexCount(n);
2026 if (numVerts == 3 || numVerts == 4) {
2028 prim.index =
Int32(n);
2030 mMesh->getIndexSpacePoint(n, 0, prim.a);
2031 mMesh->getIndexSpacePoint(n, 1, prim.b);
2032 mMesh->getIndexSpacePoint(n, 2, prim.c);
2034 evalTriangle(prim, *dataPtr);
2036 if (numVerts == 4) {
2037 mMesh->getIndexSpacePoint(n, 3, prim.b);
2038 evalTriangle(prim, *dataPtr);
2046 bool wasInterrupted()
const {
return mInterrupter && mInterrupter->wasInterrupted(); }
2048 struct Triangle { Vec3d a, b, c;
Int32 index; };
2052 enum { POLYGON_LIMIT = 1000 };
2054 SubTask(
const Triangle& prim, DataTable& dataTable,
2055 int subdivisionCount,
size_t polygonCount,
2056 Interrupter* interrupter =
nullptr)
2057 : mLocalDataTable(&dataTable)
2059 , mSubdivisionCount(subdivisionCount)
2060 , mPolygonCount(polygonCount)
2061 , mInterrupter(interrupter)
2065 void operator()()
const 2067 if (mSubdivisionCount <= 0 || mPolygonCount >= POLYGON_LIMIT) {
2069 typename VoxelizationDataType::Ptr& dataPtr = mLocalDataTable->local();
2070 if (!dataPtr) dataPtr.reset(
new VoxelizationDataType());
2072 voxelizeTriangle(mPrim, *dataPtr, mInterrupter);
2074 }
else if (!(mInterrupter && mInterrupter->wasInterrupted())) {
2075 spawnTasks(mPrim, *mLocalDataTable, mSubdivisionCount, mPolygonCount, mInterrupter);
2079 DataTable *
const mLocalDataTable;
2080 Triangle
const mPrim;
2081 int const mSubdivisionCount;
2082 size_t const mPolygonCount;
2083 Interrupter *
const mInterrupter;
2086 inline static int evalSubdivisionCount(
const Triangle& prim)
2088 const double ax = prim.a[0], bx = prim.b[0], cx = prim.c[0];
2091 const double ay = prim.a[1], by = prim.b[1], cy = prim.c[1];
2094 const double az = prim.a[2], bz = prim.b[2], cz = prim.c[2];
2097 return int(
std::max(dx,
std::max(dy, dz)) /
double(TreeType::LeafNodeType::DIM * 2));
2100 void evalTriangle(
const Triangle& prim, VoxelizationDataType& data)
const 2102 const size_t polygonCount = mMesh->polygonCount();
2103 const int subdivisionCount =
2104 polygonCount < SubTask::POLYGON_LIMIT ? evalSubdivisionCount(prim) : 0;
2106 if (subdivisionCount <= 0) {
2107 voxelizeTriangle(prim, data, mInterrupter);
2109 spawnTasks(prim, *mDataTable, subdivisionCount, polygonCount, mInterrupter);
2113 static void spawnTasks(
2114 const Triangle& mainPrim,
2115 DataTable& dataTable,
2116 int subdivisionCount,
2117 size_t polygonCount,
2118 Interrupter*
const interrupter)
2120 subdivisionCount -= 1;
2123 tbb::task_group tasks;
2125 const Vec3d ac = (mainPrim.a + mainPrim.c) * 0.5;
2126 const Vec3d bc = (mainPrim.b + mainPrim.c) * 0.5;
2127 const Vec3d ab = (mainPrim.a + mainPrim.b) * 0.5;
2130 prim.index = mainPrim.index;
2132 prim.a = mainPrim.a;
2135 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount, interrupter));
2140 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount, interrupter));
2143 prim.b = mainPrim.b;
2145 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount, interrupter));
2149 prim.c = mainPrim.c;
2150 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount, interrupter));
2155 static void voxelizeTriangle(
const Triangle& prim, VoxelizationDataType& data, Interrupter*
const interrupter)
2157 std::deque<Coord> coordList;
2160 ijk = Coord::floor(prim.a);
2161 coordList.push_back(ijk);
2166 updateDistance(ijk, prim, data);
2168 unsigned char primId = data.getNewPrimId();
2169 data.primIdAcc.setValueOnly(ijk, primId);
2171 while (!coordList.empty()) {
2172 if (interrupter && interrupter->wasInterrupted()) {
2173 thread::cancelGroupExecution();
2176 for (
Int32 pass = 0; pass < 1048576 && !coordList.empty(); ++pass) {
2177 ijk = coordList.back();
2178 coordList.pop_back();
2180 for (
Int32 i = 0; i < 26; ++i) {
2182 if (primId != data.primIdAcc.getValue(nijk)) {
2183 data.primIdAcc.setValueOnly(nijk, primId);
2184 if(updateDistance(nijk, prim, data)) coordList.push_back(nijk);
2191 static bool updateDistance(
const Coord& ijk,
const Triangle& prim, VoxelizationDataType& data)
2193 Vec3d uvw, voxelCenter(ijk[0], ijk[1], ijk[2]);
2197 const ValueType dist = ValueType((voxelCenter -
2202 if (std::isnan(dist))
2205 const ValueType oldDist = data.distAcc.getValue(ijk);
2207 if (dist < oldDist) {
2208 data.distAcc.setValue(ijk, dist);
2209 data.indexAcc.setValue(ijk, prim.index);
2213 data.indexAcc.setValueOnly(ijk,
std::min(prim.index, data.indexAcc.getValue(ijk)));
2216 return !(dist > 0.75);
2219 DataTable *
const mDataTable;
2221 Interrupter *
const mInterrupter;
2228 template<
typename TreeType>
2229 struct DiffLeafNodeMask
2234 using BoolTreeType =
typename TreeType::template ValueConverter<bool>::Type;
2235 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
2237 DiffLeafNodeMask(
const TreeType& rhsTree,
2238 std::vector<BoolLeafNodeType*>& lhsNodes)
2239 : mRhsTree(&rhsTree), mLhsNodes(lhsNodes.empty() ?
nullptr : &lhsNodes[0])
2243 void operator()(
const tbb::blocked_range<size_t>& range)
const {
2247 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2249 BoolLeafNodeType* lhsNode = mLhsNodes[n];
2250 const LeafNodeType* rhsNode = acc.
probeConstLeaf(lhsNode->origin());
2252 if (rhsNode) lhsNode->topologyDifference(*rhsNode,
false);
2258 BoolLeafNodeType **
const mLhsNodes;
2262 template<
typename LeafNodeTypeA,
typename LeafNodeTypeB>
2263 struct UnionValueMasks
2265 UnionValueMasks(std::vector<LeafNodeTypeA*>& nodesA, std::vector<LeafNodeTypeB*>& nodesB)
2266 : mNodesA(nodesA.empty() ?
nullptr : &nodesA[0])
2267 , mNodesB(nodesB.empty() ?
nullptr : &nodesB[0])
2271 void operator()(
const tbb::blocked_range<size_t>& range)
const {
2272 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2273 mNodesA[n]->topologyUnion(*mNodesB[n]);
2278 LeafNodeTypeA **
const mNodesA;
2279 LeafNodeTypeB **
const mNodesB;
2283 template<
typename TreeType>
2284 struct ConstructVoxelMask
2288 using BoolTreeType =
typename TreeType::template ValueConverter<bool>::Type;
2289 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
2291 ConstructVoxelMask(BoolTreeType& maskTree,
const TreeType& tree,
2292 std::vector<LeafNodeType*>& nodes)
2294 , mNodes(nodes.empty() ?
nullptr : &nodes[0])
2295 , mLocalMaskTree(
false)
2296 , mMaskTree(&maskTree)
2300 ConstructVoxelMask(ConstructVoxelMask& rhs,
tbb::split)
2302 , mNodes(rhs.mNodes)
2303 , mLocalMaskTree(
false)
2304 , mMaskTree(&mLocalMaskTree)
2308 void operator()(
const tbb::blocked_range<size_t>& range)
2310 using Iterator =
typename LeafNodeType::ValueOnCIter;
2315 Coord ijk, nijk, localCorod;
2318 for (
size_t n = range.begin(); n != range.end(); ++n) {
2320 LeafNodeType& node = *mNodes[n];
2322 CoordBBox bbox = node.getNodeBoundingBox();
2325 BoolLeafNodeType& maskNode = *maskAcc.
touchLeaf(node.origin());
2327 for (Iterator it = node.cbeginValueOn(); it; ++it) {
2328 ijk = it.getCoord();
2331 localCorod = LeafNodeType::offsetToLocalCoord(pos);
2333 if (localCorod[2] <
int(LeafNodeType::DIM - 1)) {
2335 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2337 nijk = ijk.offsetBy(0, 0, 1);
2341 if (localCorod[2] > 0) {
2343 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2345 nijk = ijk.offsetBy(0, 0, -1);
2349 if (localCorod[1] <
int(LeafNodeType::DIM - 1)) {
2350 npos = pos + LeafNodeType::DIM;
2351 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2353 nijk = ijk.offsetBy(0, 1, 0);
2357 if (localCorod[1] > 0) {
2358 npos = pos - LeafNodeType::DIM;
2359 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2361 nijk = ijk.offsetBy(0, -1, 0);
2365 if (localCorod[0] <
int(LeafNodeType::DIM - 1)) {
2366 npos = pos + LeafNodeType::DIM * LeafNodeType::DIM;
2367 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2369 nijk = ijk.offsetBy(1, 0, 0);
2373 if (localCorod[0] > 0) {
2374 npos = pos - LeafNodeType::DIM * LeafNodeType::DIM;
2375 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2377 nijk = ijk.offsetBy(-1, 0, 0);
2384 void join(ConstructVoxelMask& rhs) { mMaskTree->merge(*rhs.mMaskTree); }
2388 LeafNodeType **
const mNodes;
2390 BoolTreeType mLocalMaskTree;
2391 BoolTreeType *
const mMaskTree;
2396 template<
typename TreeType,
typename MeshDataAdapter>
2397 struct ExpandNarrowband
2401 using NodeMaskType =
typename LeafNodeType::NodeMaskType;
2403 using Int32LeafNodeType =
typename Int32TreeType::LeafNodeType;
2404 using BoolTreeType =
typename TreeType::template ValueConverter<bool>::Type;
2405 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
2412 Fragment() : idx(0), x(0), y(0), z(0), dist(0.0) {}
2415 : idx(idx_), x(x_), y(y_), z(z_), dist(dist_)
2419 bool operator<(
const Fragment& rhs)
const {
return idx < rhs.idx; }
2425 std::vector<BoolLeafNodeType*>& maskNodes,
2426 BoolTreeType& maskTree,
2428 Int32TreeType& indexTree,
2430 ValueType exteriorBandWidth,
2431 ValueType interiorBandWidth,
2432 ValueType voxelSize)
2433 : mMaskNodes(maskNodes.empty() ?
nullptr : &maskNodes[0])
2434 , mMaskTree(&maskTree)
2435 , mDistTree(&distTree)
2436 , mIndexTree(&indexTree)
2438 , mNewMaskTree(
false)
2440 , mUpdatedDistNodes()
2442 , mUpdatedIndexNodes()
2443 , mExteriorBandWidth(exteriorBandWidth)
2444 , mInteriorBandWidth(interiorBandWidth)
2445 , mVoxelSize(voxelSize)
2449 ExpandNarrowband(
const ExpandNarrowband& rhs,
tbb::split)
2450 : mMaskNodes(rhs.mMaskNodes)
2451 , mMaskTree(rhs.mMaskTree)
2452 , mDistTree(rhs.mDistTree)
2453 , mIndexTree(rhs.mIndexTree)
2455 , mNewMaskTree(
false)
2457 , mUpdatedDistNodes()
2459 , mUpdatedIndexNodes()
2460 , mExteriorBandWidth(rhs.mExteriorBandWidth)
2461 , mInteriorBandWidth(rhs.mInteriorBandWidth)
2462 , mVoxelSize(rhs.mVoxelSize)
2466 void join(ExpandNarrowband& rhs)
2468 mDistNodes.insert(mDistNodes.end(), rhs.mDistNodes.begin(), rhs.mDistNodes.end());
2469 mIndexNodes.insert(mIndexNodes.end(), rhs.mIndexNodes.begin(), rhs.mIndexNodes.end());
2471 mUpdatedDistNodes.insert(mUpdatedDistNodes.end(),
2472 rhs.mUpdatedDistNodes.begin(), rhs.mUpdatedDistNodes.end());
2474 mUpdatedIndexNodes.insert(mUpdatedIndexNodes.end(),
2475 rhs.mUpdatedIndexNodes.begin(), rhs.mUpdatedIndexNodes.end());
2477 mNewMaskTree.merge(rhs.mNewMaskTree);
2480 void operator()(
const tbb::blocked_range<size_t>& range)
2486 std::vector<Fragment> fragments;
2487 fragments.reserve(256);
2489 std::unique_ptr<LeafNodeType> newDistNodePt;
2490 std::unique_ptr<Int32LeafNodeType> newIndexNodePt;
2492 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2494 BoolLeafNodeType& maskNode = *mMaskNodes[n];
2495 if (maskNode.isEmpty())
continue;
2499 const Coord& origin = maskNode.origin();
2501 LeafNodeType * distNodePt = distAcc.
probeLeaf(origin);
2502 Int32LeafNodeType * indexNodePt = indexAcc.
probeLeaf(origin);
2504 assert(!distNodePt == !indexNodePt);
2506 bool usingNewNodes =
false;
2508 if (!distNodePt && !indexNodePt) {
2510 const ValueType backgroundDist = distAcc.
getValue(origin);
2512 if (!newDistNodePt.get() && !newIndexNodePt.get()) {
2513 newDistNodePt.reset(
new LeafNodeType(origin, backgroundDist));
2514 newIndexNodePt.reset(
new Int32LeafNodeType(origin, indexAcc.
getValue(origin)));
2517 if ((backgroundDist < ValueType(0.0)) !=
2518 (newDistNodePt->getValue(0) < ValueType(0.0))) {
2519 newDistNodePt->buffer().fill(backgroundDist);
2522 newDistNodePt->setOrigin(origin);
2523 newIndexNodePt->setOrigin(origin);
2526 distNodePt = newDistNodePt.get();
2527 indexNodePt = newIndexNodePt.get();
2529 usingNewNodes =
true;
2536 for (
typename BoolLeafNodeType::ValueOnIter it = maskNode.beginValueOn(); it; ++it) {
2537 bbox.expand(it.getCoord());
2542 gatherFragments(fragments, bbox, distAcc, indexAcc);
2547 bbox = maskNode.getNodeBoundingBox();
2549 bool updatedLeafNodes =
false;
2551 for (
typename BoolLeafNodeType::ValueOnIter it = maskNode.beginValueOn(); it; ++it) {
2553 const Coord ijk = it.getCoord();
2555 if (updateVoxel(ijk, 5, fragments, *distNodePt, *indexNodePt, &updatedLeafNodes)) {
2557 for (
Int32 i = 0; i < 6; ++i) {
2559 if (bbox.isInside(nijk)) {
2560 mask.setOn(BoolLeafNodeType::coordToOffset(nijk));
2566 for (
Int32 i = 6; i < 26; ++i) {
2568 if (bbox.isInside(nijk)) {
2569 mask.setOn(BoolLeafNodeType::coordToOffset(nijk));
2575 if (updatedLeafNodes) {
2578 mask -= indexNodePt->getValueMask();
2580 for (
typename NodeMaskType::OnIterator it = mask.beginOn(); it; ++it) {
2582 const Index pos = it.pos();
2583 const Coord ijk = maskNode.origin() + LeafNodeType::offsetToLocalCoord(pos);
2585 if (updateVoxel(ijk, 6, fragments, *distNodePt, *indexNodePt)) {
2586 for (
Int32 i = 0; i < 6; ++i) {
2593 if (usingNewNodes) {
2594 newDistNodePt->topologyUnion(*newIndexNodePt);
2595 mDistNodes.push_back(newDistNodePt.release());
2596 mIndexNodes.push_back(newIndexNodePt.release());
2598 mUpdatedDistNodes.push_back(distNodePt);
2599 mUpdatedIndexNodes.push_back(indexNodePt);
2607 BoolTreeType& newMaskTree() {
return mNewMaskTree; }
2609 std::vector<LeafNodeType*>& newDistNodes() {
return mDistNodes; }
2610 std::vector<LeafNodeType*>& updatedDistNodes() {
return mUpdatedDistNodes; }
2612 std::vector<Int32LeafNodeType*>& newIndexNodes() {
return mIndexNodes; }
2613 std::vector<Int32LeafNodeType*>& updatedIndexNodes() {
return mUpdatedIndexNodes; }
2619 gatherFragments(std::vector<Fragment>& fragments,
const CoordBBox& bbox,
2623 const Coord nodeMin = bbox.min() & ~(LeafNodeType::DIM - 1);
2624 const Coord nodeMax = bbox.max() & ~(LeafNodeType::DIM - 1);
2629 for (ijk[0] = nodeMin[0]; ijk[0] <= nodeMax[0]; ijk[0] += LeafNodeType::DIM) {
2630 for (ijk[1] = nodeMin[1]; ijk[1] <= nodeMax[1]; ijk[1] += LeafNodeType::DIM) {
2631 for (ijk[2] = nodeMin[2]; ijk[2] <= nodeMax[2]; ijk[2] += LeafNodeType::DIM) {
2632 if (LeafNodeType* distleaf = distAcc.
probeLeaf(ijk)) {
2635 ijk.offsetBy(LeafNodeType::DIM - 1));
2636 gatherFragments(fragments, region, *distleaf, *indexAcc.
probeLeaf(ijk));
2642 std::sort(fragments.begin(), fragments.end());
2646 gatherFragments(std::vector<Fragment>& fragments,
const CoordBBox& bbox,
2647 const LeafNodeType& distLeaf,
const Int32LeafNodeType& idxLeaf)
const 2649 const typename LeafNodeType::NodeMaskType& mask = distLeaf.getValueMask();
2650 const ValueType* distData = distLeaf.buffer().data();
2651 const Int32* idxData = idxLeaf.buffer().data();
2653 for (
int x = bbox.min()[0]; x <= bbox.max()[0]; ++x) {
2654 const Index xPos = (x & (LeafNodeType::DIM - 1u)) << (2 * LeafNodeType::LOG2DIM);
2655 for (
int y = bbox.min()[1]; y <= bbox.max()[1]; ++y) {
2656 const Index yPos = xPos + ((y & (LeafNodeType::DIM - 1u)) << LeafNodeType::LOG2DIM);
2657 for (
int z = bbox.min()[2]; z <= bbox.max()[2]; ++z) {
2658 const Index pos = yPos + (z & (LeafNodeType::DIM - 1u));
2659 if (mask.isOn(pos)) {
2660 fragments.push_back(Fragment(idxData[pos],x,y,z, std::abs(distData[pos])));
2670 computeDistance(
const Coord& ijk,
const Int32 manhattanLimit,
2671 const std::vector<Fragment>& fragments,
Int32& closestPrimIdx)
const 2673 Vec3d a, b, c, uvw, voxelCenter(ijk[0], ijk[1], ijk[2]);
2677 for (
size_t n = 0, N = fragments.size(); n < N; ++n) {
2679 const Fragment& fragment = fragments[n];
2680 if (lastIdx == fragment.idx)
continue;
2682 const Int32 dx = std::abs(fragment.x - ijk[0]);
2683 const Int32 dy = std::abs(fragment.y - ijk[1]);
2684 const Int32 dz = std::abs(fragment.z - ijk[2]);
2686 const Int32 manhattan = dx + dy + dz;
2687 if (manhattan > manhattanLimit)
continue;
2689 lastIdx = fragment.idx;
2691 const size_t polygon = size_t(lastIdx);
2693 mMesh->getIndexSpacePoint(polygon, 0, a);
2694 mMesh->getIndexSpacePoint(polygon, 1, b);
2695 mMesh->getIndexSpacePoint(polygon, 2, c);
2697 primDist = (voxelCenter -
2701 if (4 == mMesh->vertexCount(polygon)) {
2703 mMesh->getIndexSpacePoint(polygon, 3, b);
2706 a, b, c, voxelCenter, uvw)).lengthSqr();
2708 if (tmpDist < primDist) primDist = tmpDist;
2711 if (primDist < dist) {
2713 closestPrimIdx = lastIdx;
2717 return ValueType(std::sqrt(dist)) * mVoxelSize;
2723 updateVoxel(
const Coord& ijk,
const Int32 manhattanLimit,
2724 const std::vector<Fragment>& fragments,
2725 LeafNodeType& distLeaf, Int32LeafNodeType& idxLeaf,
bool* updatedLeafNodes =
nullptr)
2727 Int32 closestPrimIdx = 0;
2728 const ValueType distance = computeDistance(ijk, manhattanLimit, fragments, closestPrimIdx);
2730 const Index pos = LeafNodeType::coordToOffset(ijk);
2731 const bool inside = distLeaf.getValue(pos) < ValueType(0.0);
2733 bool activateNeighbourVoxels =
false;
2735 if (!inside && distance < mExteriorBandWidth) {
2736 if (updatedLeafNodes) *updatedLeafNodes =
true;
2737 activateNeighbourVoxels = (distance + mVoxelSize) < mExteriorBandWidth;
2738 distLeaf.setValueOnly(pos, distance);
2739 idxLeaf.setValueOn(pos, closestPrimIdx);
2740 }
else if (inside && distance < mInteriorBandWidth) {
2741 if (updatedLeafNodes) *updatedLeafNodes =
true;
2742 activateNeighbourVoxels = (distance + mVoxelSize) < mInteriorBandWidth;
2743 distLeaf.setValueOnly(pos, -distance);
2744 idxLeaf.setValueOn(pos, closestPrimIdx);
2747 return activateNeighbourVoxels;
2752 BoolLeafNodeType **
const mMaskNodes;
2753 BoolTreeType *
const mMaskTree;
2755 Int32TreeType *
const mIndexTree;
2759 BoolTreeType mNewMaskTree;
2761 std::vector<LeafNodeType*> mDistNodes, mUpdatedDistNodes;
2762 std::vector<Int32LeafNodeType*> mIndexNodes, mUpdatedIndexNodes;
2764 const ValueType mExteriorBandWidth, mInteriorBandWidth, mVoxelSize;
2768 template<
typename TreeType>
2772 AddNodes(
TreeType& tree, std::vector<LeafNodeType*>& nodes)
2773 : mTree(&tree) , mNodes(&nodes)
2777 void operator()()
const {
2779 std::vector<LeafNodeType*>& nodes = *mNodes;
2780 for (
size_t n = 0, N = nodes.size(); n < N; ++n) {
2786 std::vector<LeafNodeType*> *
const mNodes;
2790 template<
typename TreeType,
typename Int32TreeType,
typename BoolTreeType,
typename MeshDataAdapter>
2794 Int32TreeType& indexTree,
2795 BoolTreeType& maskTree,
2796 std::vector<typename BoolTreeType::LeafNodeType*>& maskNodes,
2802 ExpandNarrowband<TreeType, MeshDataAdapter> expandOp(maskNodes, maskTree,
2803 distTree, indexTree, mesh, exteriorBandWidth, interiorBandWidth, voxelSize);
2805 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, maskNodes.size()), expandOp);
2807 tbb::parallel_for(tbb::blocked_range<size_t>(0, expandOp.updatedIndexNodes().size()),
2808 UnionValueMasks<typename TreeType::LeafNodeType, typename Int32TreeType::LeafNodeType>(
2809 expandOp.updatedDistNodes(), expandOp.updatedIndexNodes()));
2811 tbb::task_group tasks;
2812 tasks.run(AddNodes<TreeType>(distTree, expandOp.newDistNodes()));
2813 tasks.run(AddNodes<Int32TreeType>(indexTree, expandOp.newIndexNodes()));
2817 maskTree.merge(expandOp.newMaskTree());
2825 template<
typename TreeType>
2826 struct TransformValues
2831 TransformValues(std::vector<LeafNodeType*>& nodes,
2832 ValueType voxelSize,
bool unsignedDist)
2834 , mVoxelSize(voxelSize)
2835 , mUnsigned(unsignedDist)
2839 void operator()(
const tbb::blocked_range<size_t>& range)
const {
2841 typename LeafNodeType::ValueOnIter iter;
2843 const bool udf = mUnsigned;
2844 const ValueType w[2] = { -mVoxelSize, mVoxelSize };
2846 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2848 for (iter = mNodes[n]->beginValueOn(); iter; ++iter) {
2849 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
2850 val = w[udf || (val < ValueType(0.0))] * std::sqrt(std::abs(val));
2856 LeafNodeType * *
const mNodes;
2857 const ValueType mVoxelSize;
2858 const bool mUnsigned;
2863 template<
typename TreeType>
2864 struct InactivateValues
2869 InactivateValues(std::vector<LeafNodeType*>& nodes,
2870 ValueType exBandWidth, ValueType inBandWidth)
2871 : mNodes(nodes.empty() ?
nullptr : &nodes[0])
2872 , mExBandWidth(exBandWidth)
2873 , mInBandWidth(inBandWidth)
2877 void operator()(
const tbb::blocked_range<size_t>& range)
const {
2879 typename LeafNodeType::ValueOnIter iter;
2880 const ValueType exVal = mExBandWidth;
2881 const ValueType inVal = -mInBandWidth;
2883 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2885 for (iter = mNodes[n]->beginValueOn(); iter; ++iter) {
2887 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
2889 const bool inside = val < ValueType(0.0);
2891 if (inside && !(val > inVal)) {
2894 }
else if (!inside && !(val < exVal)) {
2903 LeafNodeType * *
const mNodes;
2904 const ValueType mExBandWidth, mInBandWidth;
2908 template<
typename TreeType>
2914 OffsetValues(std::vector<LeafNodeType*>& nodes, ValueType offset)
2915 : mNodes(nodes.empty() ?
nullptr : &nodes[0]), mOffset(offset)
2919 void operator()(
const tbb::blocked_range<size_t>& range)
const {
2921 const ValueType offset = mOffset;
2923 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2925 typename LeafNodeType::ValueOnIter iter = mNodes[n]->beginValueOn();
2927 for (; iter; ++iter) {
2928 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
2935 LeafNodeType * *
const mNodes;
2936 const ValueType mOffset;
2940 template<
typename TreeType>
2946 Renormalize(
const TreeType& tree,
const std::vector<LeafNodeType*>& nodes,
2947 ValueType* buffer, ValueType voxelSize)
2949 , mNodes(nodes.empty() ?
nullptr : &nodes[0])
2951 , mVoxelSize(voxelSize)
2955 void operator()(
const tbb::blocked_range<size_t>& range)
const 2964 const ValueType dx = mVoxelSize, invDx = ValueType(1.0) / mVoxelSize;
2966 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2968 ValueType* bufferData = &mBuffer[n * LeafNodeType::SIZE];
2970 typename LeafNodeType::ValueOnCIter iter = mNodes[n]->cbeginValueOn();
2971 for (; iter; ++iter) {
2973 const ValueType phi0 = *iter;
2975 ijk = iter.getCoord();
2977 up[0] = acc.
getValue(ijk.offsetBy(1, 0, 0)) - phi0;
2978 up[1] = acc.
getValue(ijk.offsetBy(0, 1, 0)) - phi0;
2979 up[2] = acc.
getValue(ijk.offsetBy(0, 0, 1)) - phi0;
2981 down[0] = phi0 - acc.
getValue(ijk.offsetBy(-1, 0, 0));
2982 down[1] = phi0 - acc.
getValue(ijk.offsetBy(0, -1, 0));
2983 down[2] = phi0 - acc.
getValue(ijk.offsetBy(0, 0, -1));
2987 const ValueType diff =
math::Sqrt(normSqGradPhi) * invDx - ValueType(1.0);
2990 bufferData[iter.pos()] = phi0 - dx * S * diff;
2997 LeafNodeType
const *
const *
const mNodes;
2998 ValueType *
const mBuffer;
3000 const ValueType mVoxelSize;
3004 template<
typename TreeType>
3010 MinCombine(std::vector<LeafNodeType*>& nodes,
const ValueType* buffer)
3011 : mNodes(nodes.empty() ?
nullptr : &nodes[0]), mBuffer(buffer)
3015 void operator()(
const tbb::blocked_range<size_t>& range)
const {
3017 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
3019 const ValueType* bufferData = &mBuffer[n * LeafNodeType::SIZE];
3021 typename LeafNodeType::ValueOnIter iter = mNodes[n]->beginValueOn();
3023 for (; iter; ++iter) {
3024 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
3025 val =
std::min(val, bufferData[iter.pos()]);
3031 LeafNodeType * *
const mNodes;
3032 ValueType
const *
const mBuffer;
3046 template <
typename FloatTreeT>
3050 using ConnectivityTable = mesh_to_volume_internal::LeafNodeConnectivityTable<FloatTreeT>;
3055 ConnectivityTable nodeConnectivity(tree);
3057 std::vector<size_t> zStartNodes, yStartNodes, xStartNodes;
3062 for (
size_t n = 0; n < nodeConnectivity.size(); ++n) {
3063 if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevX()[n]) {
3064 xStartNodes.push_back(n);
3067 if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevY()[n]) {
3068 yStartNodes.push_back(n);
3071 if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevZ()[n]) {
3072 zStartNodes.push_back(n);
3076 using SweepingOp = mesh_to_volume_internal::SweepExteriorSign<FloatTreeT>;
3081 tbb::parallel_for(tbb::blocked_range<size_t>(0, zStartNodes.size()),
3084 tbb::parallel_for(tbb::blocked_range<size_t>(0, yStartNodes.size()),
3087 tbb::parallel_for(tbb::blocked_range<size_t>(0, xStartNodes.size()),
3090 const size_t numLeafNodes = nodeConnectivity.size();
3091 const size_t numVoxels = numLeafNodes * FloatTreeT::LeafNodeType::SIZE;
3093 std::unique_ptr<bool[]> changedNodeMaskA{
new bool[numLeafNodes]};
3094 std::unique_ptr<bool[]> changedNodeMaskB{
new bool[numLeafNodes]};
3095 std::unique_ptr<bool[]> changedVoxelMask{
new bool[numVoxels]};
3097 mesh_to_volume_internal::fillArray(changedNodeMaskA.get(),
true, numLeafNodes);
3098 mesh_to_volume_internal::fillArray(changedNodeMaskB.get(),
false, numLeafNodes);
3099 mesh_to_volume_internal::fillArray(changedVoxelMask.get(),
false, numVoxels);
3101 const tbb::blocked_range<size_t> nodeRange(0, numLeafNodes);
3103 bool nodesUpdated =
false;
3108 tbb::parallel_for(nodeRange, mesh_to_volume_internal::SeedFillExteriorSign<FloatTreeT>(
3109 nodeConnectivity.nodes(), changedNodeMaskA.get()));
3117 tbb::parallel_for(nodeRange, mesh_to_volume_internal::SeedPoints<FloatTreeT>(
3118 nodeConnectivity, changedNodeMaskA.get(), changedNodeMaskB.get(),
3119 changedVoxelMask.get()));
3123 changedNodeMaskA.swap(changedNodeMaskB);
3125 nodesUpdated =
false;
3126 for (
size_t n = 0; n < numLeafNodes; ++n) {
3127 nodesUpdated |= changedNodeMaskA[n];
3128 if (nodesUpdated)
break;
3134 tbb::parallel_for(nodeRange, mesh_to_volume_internal::SyncVoxelMask<FloatTreeT>(
3135 nodeConnectivity.nodes(), changedNodeMaskA.get(), changedVoxelMask.get()));
3137 }
while (nodesUpdated);
3144 template <
typename T, Index Log2Dim,
typename InteriorTest>
3167 const auto DIM = leafNode.
DIM;
3168 const auto SIZE = leafNode.
SIZE;
3170 std::vector<VoxelState> voxelState(SIZE, NOT_VISITED);
3172 std::vector<std::pair<Index, VoxelState>> offsetStack;
3173 offsetStack.reserve(SIZE);
3175 for (
Index offset=0; offset<SIZE; offset++) {
3176 const auto value = leafNode.
getValue(offset);
3180 if (std::abs(value) <= 0.75) {
3181 voxelState[offset] = NOT_ASSIGNED;
3182 }
else if (voxelState[offset] == NOT_VISITED) {
3186 if (interiorTest(coord)){
3189 offsetStack.push_back({offset, POSITIVE});
3190 voxelState[offset] = POSITIVE;
3192 offsetStack.push_back({offset, NEGATIVE});
3193 voxelState[offset] = NEGATIVE;
3196 while(!offsetStack.empty()){
3198 auto [off, state] = offsetStack[offsetStack.size()-1];
3199 offsetStack.pop_back();
3201 if (state == NEGATIVE) {
3208 for (
int dim=2; dim>=0; dim--){
3209 for (
int i = -1; i <=1; ++(++i)){
3210 int dimIdx = (off >> dim * Log2Dim) % DIM;
3211 auto neighOff = off + (1 << dim * Log2Dim) * i;
3213 (dimIdx < (int)DIM - 1) &&
3214 (voxelState[neighOff] == NOT_VISITED)) {
3216 if (std::abs(leafNode.
getValue(neighOff)) <= 0.75) {
3217 voxelState[neighOff] = NOT_ASSIGNED;
3219 offsetStack.push_back({neighOff, state});
3220 voxelState[neighOff] = state;
3250 template <
typename FloatTreeT,
typename InteriorTest>
3254 static_assert(std::is_invocable_r<bool, InteriorTest, Coord>::value,
3255 "InteriorTest has to be a function `Coord -> bool`!");
3256 static_assert(std::is_copy_constructible_v<InteriorTest>,
3257 "InteriorTest has to be copyable!");
3259 using LeafT =
typename FloatTreeT::LeafNodeType;
3263 auto op = [interiorTest](
auto& node) {
3264 using Node = std::decay_t<decltype(node)>;
3266 if constexpr (std::is_same_v<Node, LeafT>) {
3268 for (
auto iter = node.beginValueAll(); iter; ++iter) {
3269 if (!interiorTest(iter.getCoord())) {
3270 iter.setValue(-*iter);
3275 for (
auto iter = node.beginChildOff(); iter; ++iter) {
3276 if (!interiorTest(iter.getCoord())) {
3277 iter.setValue(-*iter);
3283 openvdb::tree::NodeManager nodes(tree);
3284 nodes.foreachBottomUp(op);
3289 auto op = [interiorTest](
auto& node) {
3290 using Node = std::decay_t<decltype(node)>;
3292 if constexpr (std::is_same_v<Node, LeafT>) {
3294 LeafT& leaf =
static_cast<LeafT&
>(node);
3299 for (
auto iter = node.beginChildOff(); iter; ++iter) {
3300 if (!interiorTest(iter.getCoord())) {
3301 iter.setValue(-*iter);
3307 openvdb::tree::NodeManager nodes(tree);
3308 nodes.foreachBottomUp(op);
3315 template <
typename Gr
idType,
typename MeshDataAdapter,
typename Interrupter,
typename InteriorTest>
3316 typename GridType::Ptr
3318 Interrupter& interrupter,
3321 float exteriorBandWidth,
3322 float interiorBandWidth,
3325 InteriorTest interiorTest,
3328 using GridTypePtr =
typename GridType::Ptr;
3329 using TreeType =
typename GridType::TreeType;
3330 using LeafNodeType =
typename TreeType::LeafNodeType;
3331 using ValueType =
typename GridType::ValueType;
3334 using Int32TreeType =
typename Int32GridType::TreeType;
3336 using BoolTreeType =
typename TreeType::template ValueConverter<bool>::Type;
3343 distGrid->setTransform(transform.
copy());
3345 ValueType exteriorWidth = ValueType(exteriorBandWidth);
3346 ValueType interiorWidth = ValueType(interiorBandWidth);
3350 if (!std::isfinite(exteriorWidth) || std::isnan(interiorWidth)) {
3351 std::stringstream msg;
3352 msg <<
"Illegal narrow band width: exterior = " << exteriorWidth
3353 <<
", interior = " << interiorWidth;
3358 const ValueType voxelSize = ValueType(transform.
voxelSize()[0]);
3360 if (!std::isfinite(voxelSize) ||
math::isZero(voxelSize)) {
3361 std::stringstream msg;
3362 msg <<
"Illegal transform, voxel size = " << voxelSize;
3368 exteriorWidth *= voxelSize;
3372 interiorWidth *= voxelSize;
3380 Int32GridType* indexGrid =
nullptr;
3382 typename Int32GridType::Ptr temporaryIndexGrid;
3384 if (polygonIndexGrid) {
3385 indexGrid = polygonIndexGrid;
3388 indexGrid = temporaryIndexGrid.get();
3391 indexGrid->newTree();
3392 indexGrid->setTransform(transform.
copy());
3394 if (computeSignedDistanceField) {
3398 interiorWidth = ValueType(0.0);
3401 TreeType& distTree = distGrid->tree();
3402 Int32TreeType& indexTree = indexGrid->tree();
3410 using VoxelizationDataType = mesh_to_volume_internal::VoxelizationData<TreeType>;
3411 using DataTable = tbb::enumerable_thread_specific<typename VoxelizationDataType::Ptr>;
3415 mesh_to_volume_internal::VoxelizePolygons<TreeType, MeshDataAdapter, Interrupter>;
3417 const tbb::blocked_range<size_t> polygonRange(0, mesh.polygonCount());
3419 tbb::parallel_for(polygonRange, Voxelizer(data, mesh, &interrupter));
3421 for (
typename DataTable::iterator i = data.begin(); i != data.end(); ++i) {
3422 VoxelizationDataType& dataItem = **i;
3423 mesh_to_volume_internal::combineData(
3424 distTree, indexTree, dataItem.distTree, dataItem.indexTree);
3430 if (interrupter.wasInterrupted(30))
return distGrid;
3437 if (computeSignedDistanceField) {
3440 if constexpr (std::is_same_v<InteriorTest, std::nullptr_t>) {
3442 (void) interiorTest;
3449 bool signInitializedForEveryVoxel =
3451 !std::is_same_v<InteriorTest, std::nullptr_t> &&
3455 if (!signInitializedForEveryVoxel) {
3457 std::vector<LeafNodeType*> nodes;
3458 nodes.reserve(distTree.leafCount());
3459 distTree.getNodes(nodes);
3461 const tbb::blocked_range<size_t> nodeRange(0, nodes.size());
3464 mesh_to_volume_internal::ComputeIntersectingVoxelSign<TreeType, MeshDataAdapter>;
3466 tbb::parallel_for(nodeRange, SignOp(nodes, distTree, indexTree, mesh));
3468 if (interrupter.wasInterrupted(45))
return distGrid;
3471 if (removeIntersectingVoxels) {
3473 tbb::parallel_for(nodeRange,
3474 mesh_to_volume_internal::ValidateIntersectingVoxels<TreeType>(distTree, nodes));
3476 tbb::parallel_for(nodeRange,
3477 mesh_to_volume_internal::RemoveSelfIntersectingSurface<TreeType>(
3478 nodes, distTree, indexTree));
3486 if (interrupter.wasInterrupted(50))
return distGrid;
3488 if (distTree.activeVoxelCount() == 0) {
3490 distTree.root().setBackground(exteriorWidth,
false);
3496 std::vector<LeafNodeType*> nodes;
3497 nodes.reserve(distTree.leafCount());
3498 distTree.getNodes(nodes);
3500 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3501 mesh_to_volume_internal::TransformValues<TreeType>(
3502 nodes, voxelSize, !computeSignedDistanceField));
3506 if (computeSignedDistanceField) {
3507 distTree.root().setBackground(exteriorWidth,
false);
3513 if (interrupter.wasInterrupted(54))
return distGrid;
3520 const ValueType minBandWidth = voxelSize * ValueType(2.0);
3522 if (interiorWidth > minBandWidth || exteriorWidth > minBandWidth) {
3525 BoolTreeType maskTree(
false);
3528 std::vector<LeafNodeType*> nodes;
3529 nodes.reserve(distTree.leafCount());
3530 distTree.getNodes(nodes);
3532 mesh_to_volume_internal::ConstructVoxelMask<TreeType> op(maskTree, distTree, nodes);
3533 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), op);
3539 float progress = 54.0f, step = 0.0f;
3541 2.0 * std::ceil((
std::max(interiorWidth, exteriorWidth) - minBandWidth) / voxelSize);
3543 if (estimated <
double(maxIterations)) {
3544 maxIterations = unsigned(estimated);
3545 step = 40.0f / float(maxIterations);
3548 std::vector<typename BoolTreeType::LeafNodeType*> maskNodes;
3553 if (interrupter.wasInterrupted(
int(progress)))
return distGrid;
3555 const size_t maskNodeCount = maskTree.leafCount();
3556 if (maskNodeCount == 0)
break;
3559 maskNodes.reserve(maskNodeCount);
3560 maskTree.getNodes(maskNodes);
3562 const tbb::blocked_range<size_t> range(0, maskNodes.size());
3564 tbb::parallel_for(range,
3565 mesh_to_volume_internal::DiffLeafNodeMask<TreeType>(distTree, maskNodes));
3567 mesh_to_volume_internal::expandNarrowband(distTree, indexTree, maskTree, maskNodes,
3568 mesh, exteriorWidth, interiorWidth, voxelSize);
3570 if ((++count) >= maxIterations)
break;
3575 if (interrupter.wasInterrupted(94))
return distGrid;
3577 if (!polygonIndexGrid) indexGrid->clear();
3585 if (computeSignedDistanceField && renormalizeValues) {
3587 std::vector<LeafNodeType*> nodes;
3588 nodes.reserve(distTree.leafCount());
3589 distTree.getNodes(nodes);
3591 std::unique_ptr<ValueType[]> buffer{
new ValueType[LeafNodeType::SIZE * nodes.size()]};
3593 const ValueType offset = ValueType(0.8 * voxelSize);
3595 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3596 mesh_to_volume_internal::OffsetValues<TreeType>(nodes, -offset));
3598 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3599 mesh_to_volume_internal::Renormalize<TreeType>(
3600 distTree, nodes, buffer.get(), voxelSize));
3602 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3603 mesh_to_volume_internal::MinCombine<TreeType>(nodes, buffer.get()));
3605 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3606 mesh_to_volume_internal::OffsetValues<TreeType>(
3607 nodes, offset - mesh_to_volume_internal::Tolerance<ValueType>::epsilon()));
3610 if (interrupter.wasInterrupted(99))
return distGrid;
3617 if (trimNarrowBand &&
std::min(interiorWidth, exteriorWidth) < voxelSize * ValueType(4.0)) {
3619 std::vector<LeafNodeType*> nodes;
3620 nodes.reserve(distTree.leafCount());
3621 distTree.getNodes(nodes);
3623 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3624 mesh_to_volume_internal::InactivateValues<TreeType>(
3625 nodes, exteriorWidth, computeSignedDistanceField ? interiorWidth : exteriorWidth));
3628 distTree, exteriorWidth, computeSignedDistanceField ? -interiorWidth : -exteriorWidth);
3635 template <
typename Gr
idType,
typename MeshDataAdapter,
typename InteriorTest>
3636 typename GridType::Ptr
3640 float exteriorBandWidth,
3641 float interiorBandWidth,
3648 return meshToVolume<GridType>(nullInterrupter, mesh, transform,
3649 exteriorBandWidth, interiorBandWidth, flags, polygonIndexGrid);
3660 template<
typename Gr
idType,
typename Interrupter>
3661 inline typename std::enable_if<std::is_floating_point<typename GridType::ValueType>::value,
3662 typename GridType::Ptr>::type
3664 Interrupter& interrupter,
3665 const openvdb::math::Transform& xform,
3666 const std::vector<Vec3s>& points,
3667 const std::vector<Vec3I>& triangles,
3668 const std::vector<Vec4I>& quads,
3671 bool unsignedDistanceField =
false)
3673 if (points.empty()) {
3674 return typename GridType::Ptr(
new GridType(
typename GridType::ValueType(exBandWidth)));
3677 const size_t numPoints = points.size();
3678 std::unique_ptr<Vec3s[]> indexSpacePoints{
new Vec3s[numPoints]};
3681 tbb::parallel_for(tbb::blocked_range<size_t>(0, numPoints),
3682 mesh_to_volume_internal::TransformPoints<Vec3s>(
3683 &points[0], indexSpacePoints.get(), xform));
3687 if (quads.empty()) {
3690 mesh(indexSpacePoints.get(), numPoints, &triangles[0], triangles.size());
3692 return meshToVolume<GridType>(
3693 interrupter, mesh, xform, exBandWidth, inBandWidth, conversionFlags);
3695 }
else if (triangles.empty()) {
3698 mesh(indexSpacePoints.get(), numPoints, &quads[0], quads.size());
3700 return meshToVolume<GridType>(
3701 interrupter, mesh, xform, exBandWidth, inBandWidth, conversionFlags);
3706 const size_t numPrimitives = triangles.size() + quads.size();
3707 std::unique_ptr<Vec4I[]> prims{
new Vec4I[numPrimitives]};
3709 for (
size_t n = 0, N = triangles.size(); n < N; ++n) {
3710 const Vec3I& triangle = triangles[n];
3711 Vec4I& prim = prims[n];
3712 prim[0] = triangle[0];
3713 prim[1] = triangle[1];
3714 prim[2] = triangle[2];
3718 const size_t offset = triangles.size();
3719 for (
size_t n = 0, N = quads.size(); n < N; ++n) {
3720 prims[offset + n] = quads[n];
3724 mesh(indexSpacePoints.get(), numPoints, prims.get(), numPrimitives);
3726 return meshToVolume<GridType>(interrupter, mesh, xform,
3727 exBandWidth, inBandWidth, conversionFlags);
3733 template<
typename Gr
idType,
typename Interrupter>
3734 inline typename std::enable_if<!std::is_floating_point<typename GridType::ValueType>::value,
3735 typename GridType::Ptr>::type
3739 const std::vector<Vec3s>& ,
3740 const std::vector<Vec3I>& ,
3741 const std::vector<Vec4I>& ,
3747 "mesh to volume conversion is supported only for scalar floating-point grids");
3757 template<
typename Gr
idType>
3758 typename GridType::Ptr
3760 const openvdb::math::Transform& xform,
3761 const std::vector<Vec3s>& points,
3762 const std::vector<Vec3I>& triangles,
3766 return meshToLevelSet<GridType>(nullInterrupter, xform, points, triangles, halfWidth);
3770 template<
typename Gr
idType,
typename Interrupter>
3771 typename GridType::Ptr
3773 Interrupter& interrupter,
3774 const openvdb::math::Transform& xform,
3775 const std::vector<Vec3s>& points,
3776 const std::vector<Vec3I>& triangles,
3779 std::vector<Vec4I> quads(0);
3780 return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3781 halfWidth, halfWidth);
3785 template<
typename Gr
idType>
3786 typename GridType::Ptr
3788 const openvdb::math::Transform& xform,
3789 const std::vector<Vec3s>& points,
3790 const std::vector<Vec4I>& quads,
3794 return meshToLevelSet<GridType>(nullInterrupter, xform, points, quads, halfWidth);
3798 template<
typename Gr
idType,
typename Interrupter>
3799 typename GridType::Ptr
3801 Interrupter& interrupter,
3802 const openvdb::math::Transform& xform,
3803 const std::vector<Vec3s>& points,
3804 const std::vector<Vec4I>& quads,
3807 std::vector<Vec3I> triangles(0);
3808 return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3809 halfWidth, halfWidth);
3813 template<
typename Gr
idType>
3814 typename GridType::Ptr
3816 const openvdb::math::Transform& xform,
3817 const std::vector<Vec3s>& points,
3818 const std::vector<Vec3I>& triangles,
3819 const std::vector<Vec4I>& quads,
3823 return meshToLevelSet<GridType>(
3824 nullInterrupter, xform, points, triangles, quads, halfWidth);
3828 template<
typename Gr
idType,
typename Interrupter>
3829 typename GridType::Ptr
3831 Interrupter& interrupter,
3832 const openvdb::math::Transform& xform,
3833 const std::vector<Vec3s>& points,
3834 const std::vector<Vec3I>& triangles,
3835 const std::vector<Vec4I>& quads,
3838 return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3839 halfWidth, halfWidth);
3843 template<
typename Gr
idType>
3844 typename GridType::Ptr
3846 const openvdb::math::Transform& xform,
3847 const std::vector<Vec3s>& points,
3848 const std::vector<Vec3I>& triangles,
3849 const std::vector<Vec4I>& quads,
3854 return meshToSignedDistanceField<GridType>(
3855 nullInterrupter, xform, points, triangles, quads, exBandWidth, inBandWidth);
3859 template<
typename Gr
idType,
typename Interrupter>
3860 typename GridType::Ptr
3862 Interrupter& interrupter,
3863 const openvdb::math::Transform& xform,
3864 const std::vector<Vec3s>& points,
3865 const std::vector<Vec3I>& triangles,
3866 const std::vector<Vec4I>& quads,
3870 return doMeshConversion<GridType>(interrupter, xform, points, triangles,
3871 quads, exBandWidth, inBandWidth);
3875 template<
typename Gr
idType>
3876 typename GridType::Ptr
3878 const openvdb::math::Transform& xform,
3879 const std::vector<Vec3s>& points,
3880 const std::vector<Vec3I>& triangles,
3881 const std::vector<Vec4I>& quads,
3885 return meshToUnsignedDistanceField<GridType>(
3886 nullInterrupter, xform, points, triangles, quads, bandWidth);
3890 template<
typename Gr
idType,
typename Interrupter>
3891 typename GridType::Ptr
3893 Interrupter& interrupter,
3894 const openvdb::math::Transform& xform,
3895 const std::vector<Vec3s>& points,
3896 const std::vector<Vec3I>& triangles,
3897 const std::vector<Vec4I>& quads,
3900 return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3901 bandWidth, bandWidth,
true);
3909 inline std::ostream&
3912 ostr <<
"{[ " << rhs.
mXPrim <<
", " << rhs.
mXDist <<
"]";
3913 ostr <<
" [ " << rhs.
mYPrim <<
", " << rhs.
mYDist <<
"]";
3914 ostr <<
" [ " << rhs.
mZPrim <<
", " << rhs.
mZDist <<
"]}";
3934 const std::vector<Vec3s>& pointList,
3935 const std::vector<Vec4I>& polygonList);
3937 void run(
bool threaded =
true);
3940 inline void operator() (
const tbb::blocked_range<size_t> &range);
3948 struct Primitive { Vec3d a, b, c, d;
Int32 index; };
3950 template<
bool IsQuad>
3951 inline void voxelize(
const Primitive&);
3953 template<
bool IsQuad>
3954 inline bool evalPrimitive(
const Coord&,
const Primitive&);
3956 inline bool rayTriangleIntersection(
const Vec3d& origin,
const Vec3d& dir,
3957 const Vec3d& a,
const Vec3d& b,
const Vec3d& c,
double& t);
3963 const std::vector<Vec3s>& mPointList;
3964 const std::vector<Vec4I>& mPolygonList;
3968 IntTreeT mLastPrimTree;
3974 MeshToVoxelEdgeData::GenEdgeData::GenEdgeData(
3975 const std::vector<Vec3s>& pointList,
3976 const std::vector<Vec4I>& polygonList)
3979 , mPointList(pointList)
3980 , mPolygonList(polygonList)
3982 , mLastPrimAccessor(mLastPrimTree)
3991 , mPointList(rhs.mPointList)
3992 , mPolygonList(rhs.mPolygonList)
3994 , mLastPrimAccessor(mLastPrimTree)
4003 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPolygonList.size()), *
this);
4005 (*this)(tbb::blocked_range<size_t>(0, mPolygonList.size()));
4014 using NodeChainType = RootNodeType::NodeChainType;
4015 static_assert(NodeChainType::Size > 1,
"expected tree height > 1");
4016 using InternalNodeType =
typename NodeChainType::template Get<1>;
4024 for ( ; leafIt; ++leafIt) {
4025 ijk = leafIt->origin();
4032 InternalNodeType* node = rhs.mAccessor.
getNode<InternalNodeType>();
4034 rhs.mAccessor.
clear();
4044 if (!lhsLeafPt->isValueOn(offset)) {
4045 lhsLeafPt->setValueOn(offset, rhsValue);
4077 for (
size_t n = range.begin(); n < range.end(); ++n) {
4079 const Vec4I& verts = mPolygonList[n];
4081 prim.index =
Int32(n);
4082 prim.a = Vec3d(mPointList[verts[0]]);
4083 prim.b = Vec3d(mPointList[verts[1]]);
4084 prim.c = Vec3d(mPointList[verts[2]]);
4087 prim.d = Vec3d(mPointList[verts[3]]);
4088 voxelize<true>(prim);
4090 voxelize<false>(prim);
4096 template<
bool IsQuad>
4098 MeshToVoxelEdgeData::GenEdgeData::voxelize(
const Primitive& prim)
4100 std::deque<Coord> coordList;
4103 ijk = Coord::floor(prim.a);
4104 coordList.push_back(ijk);
4106 evalPrimitive<IsQuad>(ijk, prim);
4108 while (!coordList.empty()) {
4110 ijk = coordList.back();
4111 coordList.pop_back();
4113 for (
Int32 i = 0; i < 26; ++i) {
4116 if (prim.index != mLastPrimAccessor.getValue(nijk)) {
4117 mLastPrimAccessor.setValue(nijk, prim.index);
4118 if(evalPrimitive<IsQuad>(nijk, prim)) coordList.push_back(nijk);
4125 template<
bool IsQuad>
4127 MeshToVoxelEdgeData::GenEdgeData::evalPrimitive(
const Coord& ijk,
const Primitive& prim)
4129 Vec3d uvw, org(ijk[0], ijk[1], ijk[2]);
4130 bool intersecting =
false;
4137 double dist = (org -
4140 if (rayTriangleIntersection(org, Vec3d(1.0, 0.0, 0.0), prim.a, prim.c, prim.b, t)) {
4141 if (t < edgeData.
mXDist) {
4142 edgeData.
mXDist = float(t);
4143 edgeData.
mXPrim = prim.index;
4144 intersecting =
true;
4148 if (rayTriangleIntersection(org, Vec3d(0.0, 1.0, 0.0), prim.a, prim.c, prim.b, t)) {
4149 if (t < edgeData.
mYDist) {
4150 edgeData.
mYDist = float(t);
4151 edgeData.
mYPrim = prim.index;
4152 intersecting =
true;
4156 if (rayTriangleIntersection(org, Vec3d(0.0, 0.0, 1.0), prim.a, prim.c, prim.b, t)) {
4157 if (t < edgeData.
mZDist) {
4158 edgeData.
mZDist = float(t);
4159 edgeData.
mZPrim = prim.index;
4160 intersecting =
true;
4166 double secondDist = (org -
4169 if (secondDist < dist) dist = secondDist;
4171 if (rayTriangleIntersection(org, Vec3d(1.0, 0.0, 0.0), prim.a, prim.d, prim.c, t)) {
4172 if (t < edgeData.
mXDist) {
4173 edgeData.
mXDist = float(t);
4174 edgeData.
mXPrim = prim.index;
4175 intersecting =
true;
4179 if (rayTriangleIntersection(org, Vec3d(0.0, 1.0, 0.0), prim.a, prim.d, prim.c, t)) {
4180 if (t < edgeData.
mYDist) {
4181 edgeData.
mYDist = float(t);
4182 edgeData.
mYPrim = prim.index;
4183 intersecting =
true;
4187 if (rayTriangleIntersection(org, Vec3d(0.0, 0.0, 1.0), prim.a, prim.d, prim.c, t)) {
4188 if (t < edgeData.
mZDist) {
4189 edgeData.
mZDist = float(t);
4190 edgeData.
mZPrim = prim.index;
4191 intersecting =
true;
4196 if (intersecting) mAccessor.
setValue(ijk, edgeData);
4198 return (dist < 0.86602540378443861);
4203 MeshToVoxelEdgeData::GenEdgeData::rayTriangleIntersection(
4204 const Vec3d& origin,
const Vec3d& dir,
4205 const Vec3d& a,
const Vec3d& b,
const Vec3d& c,
4212 Vec3d s1 = dir.
cross(e2);
4214 double divisor = s1.
dot(e1);
4215 if (!(std::abs(divisor) > 0.0))
return false;
4219 double inv_divisor = 1.0 / divisor;
4220 Vec3d d = origin - a;
4221 double b1 = d.
dot(s1) * inv_divisor;
4223 if (b1 < 0.0 || b1 > 1.0)
return false;
4225 Vec3d s2 = d.
cross(e1);
4226 double b2 = dir.
dot(s2) * inv_divisor;
4228 if (b2 < 0.0 || (b1 + b2) > 1.0)
return false;
4232 t = e2.dot(s2) * inv_divisor;
4233 return (t < 0.0) ?
false :
true;
4249 const std::vector<Vec3s>& pointList,
4250 const std::vector<Vec4I>& polygonList)
4264 std::vector<Vec3d>& points,
4265 std::vector<Index32>& primitives)
4275 point[0] = double(coord[0]) + data.
mXDist;
4276 point[1] = double(coord[1]);
4277 point[2] = double(coord[2]);
4279 points.push_back(point);
4280 primitives.push_back(data.
mXPrim);
4284 point[0] = double(coord[0]);
4285 point[1] = double(coord[1]) + data.
mYDist;
4286 point[2] = double(coord[2]);
4288 points.push_back(point);
4289 primitives.push_back(data.
mYPrim);
4293 point[0] = double(coord[0]);
4294 point[1] = double(coord[1]);
4295 point[2] = double(coord[2]) + data.
mZDist;
4297 points.push_back(point);
4298 primitives.push_back(data.
mZPrim);
4308 point[0] = double(coord[0]);
4309 point[1] = double(coord[1]) + data.
mYDist;
4310 point[2] = double(coord[2]);
4312 points.push_back(point);
4313 primitives.push_back(data.
mYPrim);
4317 point[0] = double(coord[0]);
4318 point[1] = double(coord[1]);
4319 point[2] = double(coord[2]) + data.
mZDist;
4321 points.push_back(point);
4322 primitives.push_back(data.
mZPrim);
4330 point[0] = double(coord[0]);
4331 point[1] = double(coord[1]) + data.
mYDist;
4332 point[2] = double(coord[2]);
4334 points.push_back(point);
4335 primitives.push_back(data.
mYPrim);
4344 point[0] = double(coord[0]) + data.
mXDist;
4345 point[1] = double(coord[1]);
4346 point[2] = double(coord[2]);
4348 points.push_back(point);
4349 primitives.push_back(data.
mXPrim);
4353 point[0] = double(coord[0]);
4354 point[1] = double(coord[1]) + data.
mYDist;
4355 point[2] = double(coord[2]);
4357 points.push_back(point);
4358 primitives.push_back(data.
mYPrim);
4368 point[0] = double(coord[0]) + data.
mXDist;
4369 point[1] = double(coord[1]);
4370 point[2] = double(coord[2]);
4372 points.push_back(point);
4373 primitives.push_back(data.
mXPrim);
4382 point[0] = double(coord[0]) + data.
mXDist;
4383 point[1] = double(coord[1]);
4384 point[2] = double(coord[2]);
4386 points.push_back(point);
4387 primitives.push_back(data.
mXPrim);
4391 point[0] = double(coord[0]);
4392 point[1] = double(coord[1]);
4393 point[2] = double(coord[2]) + data.
mZDist;
4395 points.push_back(point);
4396 primitives.push_back(data.
mZPrim);
4405 point[0] = double(coord[0]);
4406 point[1] = double(coord[1]);
4407 point[2] = double(coord[2]) + data.
mZDist;
4409 points.push_back(point);
4410 primitives.push_back(data.
mZPrim);
4416 template<
typename Gr
idType,
typename VecType>
4417 typename GridType::Ptr
4419 const openvdb::math::Transform& xform,
4420 typename VecType::ValueType halfWidth)
4426 points[0] =
Vec3s(pmin[0], pmin[1], pmin[2]);
4427 points[1] =
Vec3s(pmin[0], pmin[1], pmax[2]);
4428 points[2] =
Vec3s(pmax[0], pmin[1], pmax[2]);
4429 points[3] =
Vec3s(pmax[0], pmin[1], pmin[2]);
4430 points[4] =
Vec3s(pmin[0], pmax[1], pmin[2]);
4431 points[5] =
Vec3s(pmin[0], pmax[1], pmax[2]);
4432 points[6] =
Vec3s(pmax[0], pmax[1], pmax[2]);
4433 points[7] =
Vec3s(pmax[0], pmax[1], pmin[2]);
4436 faces[0] =
Vec4I(0, 1, 2, 3);
4437 faces[1] =
Vec4I(7, 6, 5, 4);
4438 faces[2] =
Vec4I(4, 5, 1, 0);
4439 faces[3] =
Vec4I(6, 7, 3, 2);
4440 faces[4] =
Vec4I(0, 3, 7, 4);
4441 faces[5] =
Vec4I(1, 5, 6, 2);
4445 return meshToVolume<GridType>(mesh, xform,
static_cast<float>(halfWidth), static_cast<float>(halfWidth));
4454 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION 4456 #ifdef OPENVDB_INSTANTIATE_MESHTOVOLUME 4460 #define _FUNCTION(TreeT) \ 4461 Grid<TreeT>::Ptr meshToVolume<Grid<TreeT>>(util::NullInterrupter&, \ 4462 const QuadAndTriangleDataAdapter<Vec3s, Vec3I>&, const openvdb::math::Transform&, \ 4463 float, float, int, Grid<TreeT>::ValueConverter<Int32>::Type*, std::nullptr_t, InteriorTestStrategy) 4467 #define _FUNCTION(TreeT) \ 4468 Grid<TreeT>::Ptr meshToVolume<Grid<TreeT>>(util::NullInterrupter&, \ 4469 const QuadAndTriangleDataAdapter<Vec3s, Vec4I>&, const openvdb::math::Transform&, \ 4470 float, float, int, Grid<TreeT>::ValueConverter<Int32>::Type*, std::nullptr_t, InteriorTestStrategy) 4474 #define _FUNCTION(TreeT) \ 4475 Grid<TreeT>::Ptr meshToLevelSet<Grid<TreeT>>(util::NullInterrupter&, \ 4476 const openvdb::math::Transform&, const std::vector<Vec3s>&, const std::vector<Vec3I>&, \ 4481 #define _FUNCTION(TreeT) \ 4482 Grid<TreeT>::Ptr meshToLevelSet<Grid<TreeT>>(util::NullInterrupter&, \ 4483 const openvdb::math::Transform&, const std::vector<Vec3s>&, const std::vector<Vec4I>&, \ 4488 #define _FUNCTION(TreeT) \ 4489 Grid<TreeT>::Ptr meshToLevelSet<Grid<TreeT>>(util::NullInterrupter&, \ 4490 const openvdb::math::Transform&, const std::vector<Vec3s>&, \ 4491 const std::vector<Vec3I>&, const std::vector<Vec4I>&, float) 4495 #define _FUNCTION(TreeT) \ 4496 Grid<TreeT>::Ptr meshToSignedDistanceField<Grid<TreeT>>(util::NullInterrupter&, \ 4497 const openvdb::math::Transform&, const std::vector<Vec3s>&, \ 4498 const std::vector<Vec3I>&, const std::vector<Vec4I>&, float, float) 4502 #define _FUNCTION(TreeT) \ 4503 Grid<TreeT>::Ptr meshToUnsignedDistanceField<Grid<TreeT>>(util::NullInterrupter&, \ 4504 const openvdb::math::Transform&, const std::vector<Vec3s>&, \ 4505 const std::vector<Vec3I>&, const std::vector<Vec4I>&, float) 4509 #define _FUNCTION(TreeT) \ 4510 Grid<TreeT>::Ptr createLevelSetBox<Grid<TreeT>>(const math::BBox<Vec3s>&, \ 4511 const openvdb::math::Transform&, float) 4515 #define _FUNCTION(TreeT) \ 4516 Grid<TreeT>::Ptr createLevelSetBox<Grid<TreeT>>(const math::BBox<Vec3d>&, \ 4517 const openvdb::math::Transform&, double) 4521 #define _FUNCTION(TreeT) \ 4522 void traceExteriorBoundaries(TreeT&) 4526 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION 4533 #endif // OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED Vec3< T > cross(const Vec3< T > &v) const
Return the cross product of "this" vector and v;.
Definition: Vec3.h:220
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:1543
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:49
const ValueT & getValue() const
Return the tile or voxel value to which this iterator is currently pointing.
Definition: TreeIterator.h:692
Vec2< T > minComponent(const Vec2< T > &v1, const Vec2< T > &v2)
Return component-wise minimum of the two vectors.
Definition: Vec2.h:504
The Value Accessor Implementation and API methods. The majoirty of the API matches the API of a compa...
Definition: ValueAccessor.h:68
ValueConverter<T>::Type is the type of a tree having the same hierarchy as this tree but a different ...
Definition: Tree.h:197
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
LeafNodeT * probeLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains the voxel coordinate xyz. If no LeafNode exists...
Definition: ValueAccessor.h:836
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:761
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:443
#define OPENVDB_LOG_DEBUG(message)
In debug builds only, log a debugging message of the form 'someVar << "text" << ...'.
Definition: logging.h:266
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:325
Convert polygonal meshes that consist of quads and/or triangles into signed or unsigned distance fiel...
void getNodes(ArrayT &array)
Adds all nodes of a certain type to a container with the following API:
Definition: Tree.h:1274
static const Index DIM
Definition: LeafNode.h:50
void setValue(const Coord &xyz, const ValueType &value)
Set a particular value at the given coordinate and mark the coordinate as active. ...
Definition: ValueAccessor.h:550
bool isValueOn(const Coord &xyz) const
Return the active state of the voxel at the given coordinates.
Definition: ValueAccessor.h:480
_RootNodeType RootNodeType
Definition: Tree.h:183
std::shared_ptr< T > SharedPtr
Definition: Types.h:114
Efficient multi-threaded replacement of the background values in tree.
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: LeafNode.h:1045
Base class for interrupters.
Definition: NullInterrupter.h:25
bool probeValue(const Coord &xyz, ValueType &value) const
Return the active state of the value at a given coordinate as well as its value.
Definition: ValueAccessor.h:492
static const Index SIZE
Definition: LeafNode.h:53
math::Vec4< Index32 > Vec4I
Definition: Types.h:88
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...
Base class for tree-traversal iterators over tile and voxel values.
Definition: TreeIterator.h:616
void clear()
Remove all tiles from this tree and all nodes other than the root node.
Definition: Tree.h:1297
Defined various multi-threaded utility functions for trees.
BBox< Coord > CoordBBox
Definition: NanoVDB.h:2535
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:337
void setValueOnly(const Coord &xyz, const ValueType &val)
Set the value of the voxel at the given coordinates but don't change its active state.
Definition: LeafNode.h:1103
Tree< typename RootNodeType::template ValueConverter< Int32 >::Type > Type
Definition: Tree.h:198
void clearAllAccessors()
Clear all registered accessors.
Definition: Tree.h:1356
constexpr Index32 INVALID_IDX
Definition: Util.h:19
constexpr Coord COORD_OFFSETS[26]
coordinate offset table for neighboring voxels
Definition: Util.h:22
typename RootNodeType::LeafNodeType LeafNodeType
Definition: Tree.h:186
Vec3< float > Vec3s
Definition: Vec3.h:663
Templated block class to hold specific data types and a fixed number of values determined by Log2Dim...
Definition: LeafNode.h:37
LeafIter beginLeaf()
Return an iterator over all leaf nodes in this tree.
Definition: Tree.h:975
T dot(const Vec3< T > &v) const
Dot product.
Definition: Vec3.h:191
Propagate the signs of distance values from the active voxels in the narrow band to the inactive valu...
Axis
Definition: Math.h:901
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: ValueAccessor.h:455
Definition: Exceptions.h:13
Definition: Exceptions.h:64
LeafNodeT * touchLeaf(const Coord &xyz)
Returns the leaf node that contains voxel (x, y, z) and if it doesn't exist, create it...
Definition: ValueAccessor.h:715
const Vec3T & max() const
Return a const reference to the maximum point of this bounding box.
Definition: BBox.h:64
void clear() override final
Remove all the cached nodes and invalidate the corresponding hash-keys.
Definition: ValueAccessor.h:880
Index32 Index
Definition: Types.h:54
bool operator<(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition: Tuple.h:174
void expand(ElementType padding)
Pad this bounding box.
Definition: BBox.h:321
Vec2< T > maxComponent(const Vec2< T > &v1, const Vec2< T > &v2)
Return component-wise maximum of the two vectors.
Definition: Vec2.h:513
Coord offsetToGlobalCoord(Index n) const
Return the global coordinates for a linear table offset.
Definition: LeafNode.h:1034
GridType
List of types that are currently supported by NanoVDB.
Definition: NanoVDB.h:317
const Vec3T & min() const
Return a const reference to the minimum point of this bounding box.
Definition: BBox.h:62
NodeT * getNode()
Return the node of type NodeT that has been cached on this accessor. If this accessor does not cache ...
Definition: ValueAccessor.h:848
static const Real LEVEL_SET_HALF_WIDTH
Definition: Types.h:461
OPENVDB_AX_API void run(const char *ax, openvdb::GridBase &grid, const AttributeBindings &bindings={})
Run a full AX pipeline (parse, compile and execute) on a single OpenVDB Grid.
const LeafNodeT * probeConstLeaf(const Coord &xyz) const
Return a pointer to the leaf node that contains the voxel coordinate xyz. If no LeafNode exists...
Definition: ValueAccessor.h:838
bool operator>(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition: Tuple.h:186
bool evalLeafBoundingBox(CoordBBox &bbox) const override
Return in bbox the axis-aligned bounding box of all active tiles and leaf nodes with active values...
Definition: Tree.h:1854
void setValueOn(const Coord &xyz, const ValueType &value)
Set a particular value at the given coordinate and mark the coordinate as active. ...
Definition: ValueAccessor.h:569
#define OPENVDB_REAL_TREE_INSTANTIATE(Function)
Definition: version.h.in:157
Type Pow2(Type x)
Return x2.
Definition: Math.h:548
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:729
Base class for tree-traversal iterators over all leaf nodes (but not leaf voxels) ...
Definition: TreeIterator.h:1186
RootNodeType & root()
Return this tree's root node.
Definition: Tree.h:281
Axis-aligned bounding box.
Definition: BBox.h:23
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:121
void split(ContainerT &out, const std::string &in, const char delim)
Definition: Name.h:43
typename RootNodeType::ValueType ValueType
Definition: Tree.h:184
int32_t Int32
Definition: Types.h:56
uint32_t Index32
Definition: Types.h:52
Index32 leafCount() const override
Return the number of leaf nodes.
Definition: Tree.h:340
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:212
void merge(Tree &other, MergePolicy=MERGE_ACTIVE_STATES)
Efficiently merge another tree into this tree using one of several schemes.
Definition: Tree.h:1669