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"
31 #include <tbb/blocked_range.h>
32 #include <tbb/enumerable_thread_specific.h>
33 #include <tbb/parallel_for.h>
34 #include <tbb/parallel_reduce.h>
35 #include <tbb/partitioner.h>
36 #include <tbb/task_group.h>
37 #include <tbb/task_arena.h>
45 #include <type_traits>
110 template <
typename Gr
idType,
typename MeshDataAdapter>
111 inline typename GridType::Ptr
114 const math::Transform& transform,
115 float exteriorBandWidth = 3.0f,
116 float interiorBandWidth = 3.0f,
136 template <
typename Gr
idType,
typename MeshDataAdapter,
typename Interrupter>
137 inline typename GridType::Ptr
139 Interrupter& interrupter,
141 const math::Transform& transform,
142 float exteriorBandWidth = 3.0f,
143 float interiorBandWidth = 3.0f,
161 template<
typename Po
intType,
typename PolygonType>
165 const std::vector<PolygonType>& polygons)
166 : mPointArray(points.empty() ? nullptr : &points[0])
167 , mPointArraySize(points.size())
168 , mPolygonArray(polygons.empty() ? nullptr : &polygons[0])
169 , mPolygonArraySize(polygons.size())
174 const PolygonType* polygonArray,
size_t polygonArraySize)
175 : mPointArray(pointArray)
176 , mPointArraySize(pointArraySize)
177 , mPolygonArray(polygonArray)
178 , mPolygonArraySize(polygonArraySize)
187 return (PolygonType::size == 3 || mPolygonArray[n][3] ==
util::INVALID_IDX) ? 3 : 4;
193 const PointType& p = mPointArray[mPolygonArray[n][int(v)]];
194 pos[0] = double(p[0]);
195 pos[1] = double(p[1]);
196 pos[2] = double(p[2]);
200 PointType
const *
const mPointArray;
201 size_t const mPointArraySize;
202 PolygonType
const *
const mPolygonArray;
203 size_t const mPolygonArraySize;
231 template<
typename Gr
idType>
232 inline typename GridType::Ptr
234 const openvdb::math::Transform& xform,
235 const std::vector<Vec3s>& points,
236 const std::vector<Vec3I>& triangles,
240 template<
typename Gr
idType,
typename Interrupter>
241 inline typename GridType::Ptr
243 Interrupter& interrupter,
244 const openvdb::math::Transform& xform,
245 const std::vector<Vec3s>& points,
246 const std::vector<Vec3I>& triangles,
265 template<
typename Gr
idType>
266 inline typename GridType::Ptr
268 const openvdb::math::Transform& xform,
269 const std::vector<Vec3s>& points,
270 const std::vector<Vec4I>& quads,
274 template<
typename Gr
idType,
typename Interrupter>
275 inline typename GridType::Ptr
277 Interrupter& interrupter,
278 const openvdb::math::Transform& xform,
279 const std::vector<Vec3s>& points,
280 const std::vector<Vec4I>& quads,
300 template<
typename Gr
idType>
301 inline typename GridType::Ptr
303 const openvdb::math::Transform& xform,
304 const std::vector<Vec3s>& points,
305 const std::vector<Vec3I>& triangles,
306 const std::vector<Vec4I>& quads,
310 template<
typename Gr
idType,
typename Interrupter>
311 inline typename GridType::Ptr
313 Interrupter& interrupter,
314 const openvdb::math::Transform& xform,
315 const std::vector<Vec3s>& points,
316 const std::vector<Vec3I>& triangles,
317 const std::vector<Vec4I>& quads,
339 template<
typename Gr
idType>
340 inline typename GridType::Ptr
342 const openvdb::math::Transform& xform,
343 const std::vector<Vec3s>& points,
344 const std::vector<Vec3I>& triangles,
345 const std::vector<Vec4I>& quads,
350 template<
typename Gr
idType,
typename Interrupter>
351 inline typename GridType::Ptr
353 Interrupter& interrupter,
354 const openvdb::math::Transform& xform,
355 const std::vector<Vec3s>& points,
356 const std::vector<Vec3I>& triangles,
357 const std::vector<Vec4I>& quads,
376 template<
typename Gr
idType>
377 inline typename GridType::Ptr
379 const openvdb::math::Transform& xform,
380 const std::vector<Vec3s>& points,
381 const std::vector<Vec3I>& triangles,
382 const std::vector<Vec4I>& quads,
386 template<
typename Gr
idType,
typename Interrupter>
387 inline typename GridType::Ptr
389 Interrupter& interrupter,
390 const openvdb::math::Transform& xform,
391 const std::vector<Vec3s>& points,
392 const std::vector<Vec3I>& triangles,
393 const std::vector<Vec4I>& quads,
406 template<
typename Gr
idType,
typename VecType>
407 inline typename GridType::Ptr
409 const openvdb::math::Transform& xform,
422 template <
typename FloatTreeT>
440 : mXDist(dist), mYDist(dist), mZDist(dist)
483 void convert(
const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList);
489 std::vector<Vec3d>& points, std::vector<Index32>& primitives);
510 namespace mesh_to_volume_internal {
512 template<
typename Po
intType>
513 struct TransformPoints {
515 TransformPoints(
const PointType* pointsIn, PointType* pointsOut,
517 : mPointsIn(pointsIn), mPointsOut(pointsOut), mXform(&xform)
521 void operator()(
const tbb::blocked_range<size_t>& range)
const {
525 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
527 const PointType& wsP = mPointsIn[n];
528 pos[0] = double(wsP[0]);
529 pos[1] = double(wsP[1]);
530 pos[2] = double(wsP[2]);
532 pos = mXform->worldToIndex(pos);
534 PointType& isP = mPointsOut[n];
535 isP[0] =
typename PointType::value_type(pos[0]);
536 isP[1] =
typename PointType::value_type(pos[1]);
537 isP[2] =
typename PointType::value_type(pos[2]);
541 PointType
const *
const mPointsIn;
542 PointType *
const mPointsOut;
543 math::Transform
const *
const mXform;
547 template<
typename ValueType>
550 static ValueType epsilon() {
return ValueType(1e-7); }
551 static ValueType minNarrowBandWidth() {
return ValueType(1.0 + 1e-6); }
558 template<
typename TreeType>
559 class CombineLeafNodes
565 using LeafNodeType =
typename TreeType::LeafNodeType;
566 using Int32LeafNodeType =
typename Int32TreeType::LeafNodeType;
568 CombineLeafNodes(TreeType& lhsDistTree, Int32TreeType& lhsIdxTree,
569 LeafNodeType ** rhsDistNodes, Int32LeafNodeType ** rhsIdxNodes)
570 : mDistTree(&lhsDistTree)
571 , mIdxTree(&lhsIdxTree)
572 , mRhsDistNodes(rhsDistNodes)
573 , mRhsIdxNodes(rhsIdxNodes)
577 void operator()(
const tbb::blocked_range<size_t>& range)
const {
579 tree::ValueAccessor<TreeType> distAcc(*mDistTree);
580 tree::ValueAccessor<Int32TreeType> idxAcc(*mIdxTree);
582 using DistValueType =
typename LeafNodeType::ValueType;
583 using IndexValueType =
typename Int32LeafNodeType::ValueType;
585 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
587 const Coord& origin = mRhsDistNodes[n]->origin();
589 LeafNodeType* lhsDistNode = distAcc.probeLeaf(origin);
590 Int32LeafNodeType* lhsIdxNode = idxAcc.probeLeaf(origin);
592 DistValueType* lhsDistData = lhsDistNode->buffer().data();
593 IndexValueType* lhsIdxData = lhsIdxNode->buffer().data();
595 const DistValueType* rhsDistData = mRhsDistNodes[n]->buffer().data();
596 const IndexValueType* rhsIdxData = mRhsIdxNodes[n]->buffer().data();
599 for (
Index32 offset = 0; offset < LeafNodeType::SIZE; ++offset) {
603 const DistValueType& lhsValue = lhsDistData[offset];
604 const DistValueType& rhsValue = rhsDistData[offset];
606 if (rhsValue < lhsValue) {
607 lhsDistNode->setValueOn(offset, rhsValue);
608 lhsIdxNode->setValueOn(offset, rhsIdxData[offset]);
610 lhsIdxNode->setValueOn(offset,
611 std::min(lhsIdxData[offset], rhsIdxData[offset]));
616 delete mRhsDistNodes[n];
617 delete mRhsIdxNodes[n];
623 TreeType *
const mDistTree;
624 Int32TreeType *
const mIdxTree;
626 LeafNodeType **
const mRhsDistNodes;
627 Int32LeafNodeType **
const mRhsIdxNodes;
634 template<
typename TreeType>
635 struct StashOriginAndStoreOffset
637 using LeafNodeType =
typename TreeType::LeafNodeType;
639 StashOriginAndStoreOffset(std::vector<LeafNodeType*>& nodes, Coord* coordinates)
640 : mNodes(nodes.empty() ? nullptr : &nodes[0]), mCoordinates(coordinates)
644 void operator()(
const tbb::blocked_range<size_t>& range)
const {
645 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
646 Coord& origin =
const_cast<Coord&
>(mNodes[n]->origin());
647 mCoordinates[n] = origin;
648 origin[0] =
static_cast<int>(n);
652 LeafNodeType **
const mNodes;
653 Coord *
const mCoordinates;
657 template<
typename TreeType>
660 using LeafNodeType =
typename TreeType::LeafNodeType;
662 RestoreOrigin(std::vector<LeafNodeType*>& nodes,
const Coord* coordinates)
663 : mNodes(nodes.empty() ? nullptr : &nodes[0]), mCoordinates(coordinates)
667 void operator()(
const tbb::blocked_range<size_t>& range)
const {
668 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
669 Coord& origin =
const_cast<Coord&
>(mNodes[n]->origin());
670 origin[0] = mCoordinates[n][0];
674 LeafNodeType **
const mNodes;
675 Coord
const *
const mCoordinates;
679 template<
typename TreeType>
680 class ComputeNodeConnectivity
683 using LeafNodeType =
typename TreeType::LeafNodeType;
685 ComputeNodeConnectivity(
const TreeType& tree,
const Coord* coordinates,
686 size_t* offsets,
size_t numNodes,
const CoordBBox& bbox)
688 , mCoordinates(coordinates)
690 , mNumNodes(numNodes)
695 ComputeNodeConnectivity(
const ComputeNodeConnectivity&) =
default;
698 ComputeNodeConnectivity& operator=(
const ComputeNodeConnectivity&) =
delete;
700 void operator()(
const tbb::blocked_range<size_t>& range)
const {
702 size_t* offsetsNextX = mOffsets;
703 size_t* offsetsPrevX = mOffsets + mNumNodes;
704 size_t* offsetsNextY = mOffsets + mNumNodes * 2;
705 size_t* offsetsPrevY = mOffsets + mNumNodes * 3;
706 size_t* offsetsNextZ = mOffsets + mNumNodes * 4;
707 size_t* offsetsPrevZ = mOffsets + mNumNodes * 5;
709 tree::ValueAccessor<const TreeType> acc(*mTree);
711 const Int32 DIM =
static_cast<Int32>(LeafNodeType::DIM);
713 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
714 const Coord& origin = mCoordinates[n];
715 offsetsNextX[n] = findNeighbourNode(acc, origin, Coord(DIM, 0, 0));
716 offsetsPrevX[n] = findNeighbourNode(acc, origin, Coord(-DIM, 0, 0));
717 offsetsNextY[n] = findNeighbourNode(acc, origin, Coord(0, DIM, 0));
718 offsetsPrevY[n] = findNeighbourNode(acc, origin, Coord(0, -DIM, 0));
719 offsetsNextZ[n] = findNeighbourNode(acc, origin, Coord(0, 0, DIM));
720 offsetsPrevZ[n] = findNeighbourNode(acc, origin, Coord(0, 0, -DIM));
724 size_t findNeighbourNode(tree::ValueAccessor<const TreeType>& acc,
725 const Coord& start,
const Coord& step)
const
727 Coord ijk = start + step;
728 CoordBBox bbox(mBBox);
730 while (bbox.isInside(ijk)) {
731 const LeafNodeType* node = acc.probeConstLeaf(ijk);
732 if (node)
return static_cast<size_t>(node->origin()[0]);
741 TreeType
const *
const mTree;
742 Coord
const *
const mCoordinates;
743 size_t *
const mOffsets;
745 const size_t mNumNodes;
746 const CoordBBox mBBox;
750 template<
typename TreeType>
751 struct LeafNodeConnectivityTable
755 using LeafNodeType =
typename TreeType::LeafNodeType;
757 LeafNodeConnectivityTable(TreeType& tree)
759 mLeafNodes.reserve(tree.leafCount());
760 tree.getNodes(mLeafNodes);
762 if (mLeafNodes.empty())
return;
765 tree.evalLeafBoundingBox(bbox);
767 const tbb::blocked_range<size_t> range(0, mLeafNodes.size());
771 std::unique_ptr<Coord[]> coordinates{
new Coord[mLeafNodes.size()]};
772 tbb::parallel_for(range,
773 StashOriginAndStoreOffset<TreeType>(mLeafNodes, coordinates.get()));
776 mOffsets.reset(
new size_t[mLeafNodes.size() * 6]);
779 tbb::parallel_for(range, ComputeNodeConnectivity<TreeType>(
780 tree, coordinates.get(), mOffsets.get(), mLeafNodes.size(), bbox));
783 tbb::parallel_for(range, RestoreOrigin<TreeType>(mLeafNodes, coordinates.get()));
786 size_t size()
const {
return mLeafNodes.size(); }
788 std::vector<LeafNodeType*>& nodes() {
return mLeafNodes; }
789 const std::vector<LeafNodeType*>& nodes()
const {
return mLeafNodes; }
792 const size_t* offsetsNextX()
const {
return mOffsets.get(); }
793 const size_t* offsetsPrevX()
const {
return mOffsets.get() + mLeafNodes.size(); }
795 const size_t* offsetsNextY()
const {
return mOffsets.get() + mLeafNodes.size() * 2; }
796 const size_t* offsetsPrevY()
const {
return mOffsets.get() + mLeafNodes.size() * 3; }
798 const size_t* offsetsNextZ()
const {
return mOffsets.get() + mLeafNodes.size() * 4; }
799 const size_t* offsetsPrevZ()
const {
return mOffsets.get() + mLeafNodes.size() * 5; }
802 std::vector<LeafNodeType*> mLeafNodes;
803 std::unique_ptr<size_t[]> mOffsets;
807 template<
typename TreeType>
808 class SweepExteriorSign
814 using ValueType =
typename TreeType::ValueType;
815 using LeafNodeType =
typename TreeType::LeafNodeType;
816 using ConnectivityTable = LeafNodeConnectivityTable<TreeType>;
818 SweepExteriorSign(
Axis axis,
const std::vector<size_t>& startNodeIndices,
819 ConnectivityTable& connectivity)
820 : mStartNodeIndices(startNodeIndices.empty() ? nullptr : &startNodeIndices[0])
821 , mConnectivity(&connectivity)
826 void operator()(
const tbb::blocked_range<size_t>& range)
const {
828 constexpr
Int32 DIM =
static_cast<Int32>(LeafNodeType::DIM);
830 std::vector<LeafNodeType*>& nodes = mConnectivity->nodes();
833 size_t idxA = 0, idxB = 1;
836 const size_t* nextOffsets = mConnectivity->offsetsNextZ();
837 const size_t* prevOffsets = mConnectivity->offsetsPrevZ();
845 nextOffsets = mConnectivity->offsetsNextY();
846 prevOffsets = mConnectivity->offsetsPrevY();
848 }
else if (mAxis ==
X_AXIS) {
854 nextOffsets = mConnectivity->offsetsNextX();
855 prevOffsets = mConnectivity->offsetsPrevX();
860 Int32& a = ijk[idxA];
861 Int32& b = ijk[idxB];
863 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
865 size_t startOffset = mStartNodeIndices[n];
866 size_t lastOffset = startOffset;
870 for (a = 0; a < DIM; ++a) {
871 for (b = 0; b < DIM; ++b) {
873 pos =
static_cast<Int32>(LeafNodeType::coordToOffset(ijk));
874 size_t offset = startOffset;
877 while ( offset != ConnectivityTable::INVALID_OFFSET &&
878 traceVoxelLine(*nodes[offset], pos, step) ) {
881 offset = nextOffsets[offset];
886 while (offset != ConnectivityTable::INVALID_OFFSET) {
888 offset = nextOffsets[offset];
893 pos += step * (DIM - 1);
894 while ( offset != ConnectivityTable::INVALID_OFFSET &&
895 traceVoxelLine(*nodes[offset], pos, -step)) {
896 offset = prevOffsets[offset];
904 bool traceVoxelLine(LeafNodeType& node,
Int32 pos,
const Int32 step)
const {
906 ValueType* data = node.buffer().data();
908 bool isOutside =
true;
910 for (
Index i = 0; i < LeafNodeType::DIM; ++i) {
913 ValueType& dist = data[pos];
915 if (dist < ValueType(0.0)) {
919 if (!(dist > ValueType(0.75))) isOutside =
false;
921 if (isOutside) dist = ValueType(-dist);
932 size_t const *
const mStartNodeIndices;
933 ConnectivityTable *
const mConnectivity;
939 template<
typename LeafNodeType>
941 seedFill(LeafNodeType& node)
943 using ValueType =
typename LeafNodeType::ValueType;
944 using Queue = std::deque<Index>;
947 ValueType* data = node.buffer().data();
951 for (
Index pos = 0; pos < LeafNodeType::SIZE; ++pos) {
952 if (data[pos] < 0.0) seedPoints.push_back(pos);
955 if (seedPoints.empty())
return;
958 for (Queue::iterator it = seedPoints.begin(); it != seedPoints.end(); ++it) {
959 ValueType& dist = data[*it];
966 Index pos(0), nextPos(0);
968 while (!seedPoints.empty()) {
970 pos = seedPoints.back();
971 seedPoints.pop_back();
973 ValueType& dist = data[pos];
975 if (!(dist < ValueType(0.0))) {
979 ijk = LeafNodeType::offsetToLocalCoord(pos);
982 nextPos = pos - LeafNodeType::DIM * LeafNodeType::DIM;
983 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
986 if (ijk[0] != (LeafNodeType::DIM - 1)) {
987 nextPos = pos + LeafNodeType::DIM * LeafNodeType::DIM;
988 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
992 nextPos = pos - LeafNodeType::DIM;
993 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
996 if (ijk[1] != (LeafNodeType::DIM - 1)) {
997 nextPos = pos + LeafNodeType::DIM;
998 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1003 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1006 if (ijk[2] != (LeafNodeType::DIM - 1)) {
1008 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1015 template<
typename LeafNodeType>
1017 scanFill(LeafNodeType& node)
1019 bool updatedNode =
false;
1021 using ValueType =
typename LeafNodeType::ValueType;
1022 ValueType* data = node.buffer().data();
1026 bool updatedSign =
true;
1027 while (updatedSign) {
1029 updatedSign =
false;
1031 for (
Index pos = 0; pos < LeafNodeType::SIZE; ++pos) {
1033 ValueType& dist = data[pos];
1035 if (!(dist < ValueType(0.0)) && dist > ValueType(0.75)) {
1037 ijk = LeafNodeType::offsetToLocalCoord(pos);
1040 if (ijk[2] != 0 && data[pos - 1] < ValueType(0.0)) {
1042 dist = ValueType(-dist);
1045 }
else if (ijk[2] != (LeafNodeType::DIM - 1) && data[pos + 1] < ValueType(0.0)) {
1047 dist = ValueType(-dist);
1050 }
else if (ijk[1] != 0 && data[pos - LeafNodeType::DIM] < ValueType(0.0)) {
1052 dist = ValueType(-dist);
1055 }
else if (ijk[1] != (LeafNodeType::DIM - 1)
1056 && data[pos + LeafNodeType::DIM] < ValueType(0.0))
1059 dist = ValueType(-dist);
1062 }
else if (ijk[0] != 0
1063 && data[pos - LeafNodeType::DIM * LeafNodeType::DIM] < ValueType(0.0))
1066 dist = ValueType(-dist);
1069 }
else if (ijk[0] != (LeafNodeType::DIM - 1)
1070 && data[pos + LeafNodeType::DIM * LeafNodeType::DIM] < ValueType(0.0))
1073 dist = ValueType(-dist);
1078 updatedNode |= updatedSign;
1085 template<
typename TreeType>
1086 class SeedFillExteriorSign
1089 using ValueType =
typename TreeType::ValueType;
1090 using LeafNodeType =
typename TreeType::LeafNodeType;
1092 SeedFillExteriorSign(std::vector<LeafNodeType*>& nodes,
const bool* changedNodeMask)
1093 : mNodes(nodes.empty() ? nullptr : &nodes[0])
1094 , mChangedNodeMask(changedNodeMask)
1098 void operator()(
const tbb::blocked_range<size_t>& range)
const {
1099 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1100 if (mChangedNodeMask[n]) {
1106 scanFill(*mNodes[n]);
1111 LeafNodeType **
const mNodes;
1112 const bool *
const mChangedNodeMask;
1116 template<
typename ValueType>
1119 FillArray(ValueType* array,
const ValueType v) : mArray(array), mValue(v) { }
1121 void operator()(
const tbb::blocked_range<size_t>& range)
const {
1122 const ValueType v = mValue;
1123 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1128 ValueType *
const mArray;
1129 const ValueType mValue;
1133 template<
typename ValueType>
1135 fillArray(ValueType* array,
const ValueType val,
const size_t length)
1137 const auto grainSize = std::max<size_t>(
1138 length / tbb::this_task_arena::max_concurrency(), 1024);
1139 const tbb::blocked_range<size_t> range(0, length, grainSize);
1140 tbb::parallel_for(range, FillArray<ValueType>(array, val), tbb::simple_partitioner());
1144 template<
typename TreeType>
1148 using ValueType =
typename TreeType::ValueType;
1149 using LeafNodeType =
typename TreeType::LeafNodeType;
1151 SyncVoxelMask(std::vector<LeafNodeType*>& nodes,
1152 const bool* changedNodeMask,
bool* changedVoxelMask)
1153 : mNodes(nodes.empty() ? nullptr : &nodes[0])
1154 , mChangedNodeMask(changedNodeMask)
1155 , mChangedVoxelMask(changedVoxelMask)
1159 void operator()(
const tbb::blocked_range<size_t>& range)
const {
1160 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1162 if (mChangedNodeMask[n]) {
1163 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1165 ValueType* data = mNodes[n]->buffer().data();
1167 for (
Index pos = 0; pos < LeafNodeType::SIZE; ++pos) {
1169 data[pos] = ValueType(-data[pos]);
1177 LeafNodeType **
const mNodes;
1178 bool const *
const mChangedNodeMask;
1179 bool *
const mChangedVoxelMask;
1183 template<
typename TreeType>
1187 using ValueType =
typename TreeType::ValueType;
1188 using LeafNodeType =
typename TreeType::LeafNodeType;
1189 using ConnectivityTable = LeafNodeConnectivityTable<TreeType>;
1191 SeedPoints(ConnectivityTable& connectivity,
1192 bool* changedNodeMask,
bool* nodeMask,
bool* changedVoxelMask)
1193 : mConnectivity(&connectivity)
1194 , mChangedNodeMask(changedNodeMask)
1195 , mNodeMask(nodeMask)
1196 , mChangedVoxelMask(changedVoxelMask)
1200 void operator()(
const tbb::blocked_range<size_t>& range)
const {
1202 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1203 bool changedValue =
false;
1205 changedValue |= processZ(n,
true);
1206 changedValue |= processZ(n,
false);
1208 changedValue |= processY(n,
true);
1209 changedValue |= processY(n,
false);
1211 changedValue |= processX(n,
true);
1212 changedValue |= processX(n,
false);
1214 mNodeMask[n] = changedValue;
1219 bool processZ(
const size_t n,
bool firstFace)
const
1221 const size_t offset =
1222 firstFace ? mConnectivity->offsetsPrevZ()[n] : mConnectivity->offsetsNextZ()[n];
1223 if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1225 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1227 const ValueType* lhsData = mConnectivity->nodes()[n]->buffer().data();
1228 const ValueType* rhsData = mConnectivity->nodes()[offset]->buffer().data();
1230 const Index lastOffset = LeafNodeType::DIM - 1;
1231 const Index lhsOffset =
1232 firstFace ? 0 : lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1234 Index tmpPos(0), pos(0);
1235 bool changedValue =
false;
1237 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
1238 tmpPos = x << (2 * LeafNodeType::LOG2DIM);
1239 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
1240 pos = tmpPos + (y << LeafNodeType::LOG2DIM);
1242 if (lhsData[pos + lhsOffset] > ValueType(0.75)) {
1243 if (rhsData[pos + rhsOffset] < ValueType(0.0)) {
1244 changedValue =
true;
1245 mask[pos + lhsOffset] =
true;
1251 return changedValue;
1257 bool processY(
const size_t n,
bool firstFace)
const
1259 const size_t offset =
1260 firstFace ? mConnectivity->offsetsPrevY()[n] : mConnectivity->offsetsNextY()[n];
1261 if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1263 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1265 const ValueType* lhsData = mConnectivity->nodes()[n]->buffer().data();
1266 const ValueType* rhsData = mConnectivity->nodes()[offset]->buffer().data();
1268 const Index lastOffset = LeafNodeType::DIM * (LeafNodeType::DIM - 1);
1269 const Index lhsOffset =
1270 firstFace ? 0 : lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1272 Index tmpPos(0), pos(0);
1273 bool changedValue =
false;
1275 for (
Index x = 0; x < LeafNodeType::DIM; ++x) {
1276 tmpPos = x << (2 * LeafNodeType::LOG2DIM);
1277 for (
Index z = 0; z < LeafNodeType::DIM; ++z) {
1280 if (lhsData[pos + lhsOffset] > ValueType(0.75)) {
1281 if (rhsData[pos + rhsOffset] < ValueType(0.0)) {
1282 changedValue =
true;
1283 mask[pos + lhsOffset] =
true;
1289 return changedValue;
1295 bool processX(
const size_t n,
bool firstFace)
const
1297 const size_t offset =
1298 firstFace ? mConnectivity->offsetsPrevX()[n] : mConnectivity->offsetsNextX()[n];
1299 if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1301 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1303 const ValueType* lhsData = mConnectivity->nodes()[n]->buffer().data();
1304 const ValueType* rhsData = mConnectivity->nodes()[offset]->buffer().data();
1306 const Index lastOffset = LeafNodeType::DIM * LeafNodeType::DIM * (LeafNodeType::DIM-1);
1307 const Index lhsOffset =
1308 firstFace ? 0 : lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1310 Index tmpPos(0), pos(0);
1311 bool changedValue =
false;
1313 for (
Index y = 0; y < LeafNodeType::DIM; ++y) {
1314 tmpPos = y << LeafNodeType::LOG2DIM;
1315 for (
Index z = 0; z < LeafNodeType::DIM; ++z) {
1318 if (lhsData[pos + lhsOffset] > ValueType(0.75)) {
1319 if (rhsData[pos + rhsOffset] < ValueType(0.0)) {
1320 changedValue =
true;
1321 mask[pos + lhsOffset] =
true;
1327 return changedValue;
1333 ConnectivityTable *
const mConnectivity;
1334 bool *
const mChangedNodeMask;
1335 bool *
const mNodeMask;
1336 bool *
const mChangedVoxelMask;
1342 template<
typename TreeType,
typename MeshDataAdapter>
1343 struct ComputeIntersectingVoxelSign
1345 using ValueType =
typename TreeType::ValueType;
1346 using LeafNodeType =
typename TreeType::LeafNodeType;
1348 using Int32LeafNodeType =
typename Int32TreeType::LeafNodeType;
1351 using MaskArray = std::unique_ptr<bool[]>;
1352 using LocalData = std::pair<PointArray, MaskArray>;
1353 using LocalDataTable = tbb::enumerable_thread_specific<LocalData>;
1355 ComputeIntersectingVoxelSign(
1356 std::vector<LeafNodeType*>& distNodes,
1357 const TreeType& distTree,
1358 const Int32TreeType& indexTree,
1360 : mDistNodes(distNodes.empty() ? nullptr : &distNodes[0])
1361 , mDistTree(&distTree)
1362 , mIndexTree(&indexTree)
1364 , mLocalDataTable(new LocalDataTable())
1369 void operator()(
const tbb::blocked_range<size_t>& range)
const {
1371 tree::ValueAccessor<const TreeType> distAcc(*mDistTree);
1372 tree::ValueAccessor<const Int32TreeType> idxAcc(*mIndexTree);
1376 Index xPos(0), yPos(0);
1377 Coord ijk, nijk, nodeMin, nodeMax;
1378 Vec3d cp, xyz, nxyz, dir1, dir2;
1380 LocalData& localData = mLocalDataTable->local();
1383 if (!points) points.reset(
new Vec3d[LeafNodeType::SIZE * 2]);
1385 MaskArray& mask = localData.second;
1386 if (!mask) mask.reset(
new bool[LeafNodeType::SIZE]);
1389 typename LeafNodeType::ValueOnCIter it;
1391 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1393 LeafNodeType& node = *mDistNodes[n];
1394 ValueType* data = node.buffer().data();
1396 const Int32LeafNodeType* idxNode = idxAcc.probeConstLeaf(node.origin());
1397 const Int32* idxData = idxNode->buffer().data();
1399 nodeMin = node.origin();
1400 nodeMax = nodeMin.offsetBy(LeafNodeType::DIM - 1);
1403 memset(mask.get(), 0,
sizeof(
bool) * LeafNodeType::SIZE);
1405 for (it = node.cbeginValueOn(); it; ++it) {
1406 Index pos = it.pos();
1408 ValueType& dist = data[pos];
1409 if (dist < 0.0 || dist > 0.75)
continue;
1411 ijk = node.offsetToGlobalCoord(pos);
1413 xyz[0] = double(ijk[0]);
1414 xyz[1] = double(ijk[1]);
1415 xyz[2] = double(ijk[2]);
1421 bool flipSign =
false;
1423 for (nijk[0] = bbox.min()[0]; nijk[0] <= bbox.max()[0] && !flipSign; ++nijk[0]) {
1424 xPos = (nijk[0] & (LeafNodeType::DIM - 1u)) << (2 * LeafNodeType::LOG2DIM);
1425 for (nijk[1]=bbox.min()[1]; nijk[1] <= bbox.max()[1] && !flipSign; ++nijk[1]) {
1426 yPos = xPos + ((nijk[1] & (LeafNodeType::DIM-1u)) << LeafNodeType::LOG2DIM);
1427 for (nijk[2] = bbox.min()[2]; nijk[2] <= bbox.max()[2]; ++nijk[2]) {
1428 pos = yPos + (nijk[2] & (LeafNodeType::DIM - 1u));
1430 const Int32& polyIdx = idxData[pos];
1435 const Index pointIndex = pos * 2;
1441 nxyz[0] = double(nijk[0]);
1442 nxyz[1] = double(nijk[1]);
1443 nxyz[2] = double(nijk[2]);
1445 Vec3d& point = points[pointIndex];
1447 point = closestPoint(nxyz, polyIdx);
1449 Vec3d& direction = points[pointIndex + 1];
1450 direction = nxyz - point;
1451 direction.normalize();
1454 dir1 = xyz - points[pointIndex];
1457 if (points[pointIndex + 1].dot(dir1) > 0.0) {
1468 for (
Int32 m = 0; m < 26; ++m) {
1471 if (!bbox.isInside(nijk) && distAcc.probeValue(nijk, nval) && nval<-0.75) {
1472 nxyz[0] = double(nijk[0]);
1473 nxyz[1] = double(nijk[1]);
1474 nxyz[2] = double(nijk[2]);
1476 cp = closestPoint(nxyz, idxAcc.getValue(nijk));
1484 if (dir2.dot(dir1) > 0.0) {
1500 Vec3d a, b, c, cp, uvw;
1502 const size_t polygon = size_t(polyIdx);
1503 mMesh->getIndexSpacePoint(polygon, 0, a);
1504 mMesh->getIndexSpacePoint(polygon, 1, b);
1505 mMesh->getIndexSpacePoint(polygon, 2, c);
1509 if (4 == mMesh->vertexCount(polygon)) {
1511 mMesh->getIndexSpacePoint(polygon, 3, b);
1515 if ((center - c).lengthSqr() < (center - cp).lengthSqr()) {
1524 LeafNodeType **
const mDistNodes;
1525 TreeType
const *
const mDistTree;
1526 Int32TreeType
const *
const mIndexTree;
1529 SharedPtr<LocalDataTable> mLocalDataTable;
1536 template<
typename LeafNodeType>
1538 maskNodeInternalNeighbours(
const Index pos,
bool (&mask)[26])
1540 using NodeT = LeafNodeType;
1542 const Coord ijk = NodeT::offsetToLocalCoord(pos);
1546 mask[0] = ijk[0] != (NodeT::DIM - 1);
1548 mask[1] = ijk[0] != 0;
1550 mask[2] = ijk[1] != (NodeT::DIM - 1);
1552 mask[3] = ijk[1] != 0;
1554 mask[4] = ijk[2] != (NodeT::DIM - 1);
1556 mask[5] = ijk[2] != 0;
1560 mask[6] = mask[0] && mask[5];
1562 mask[7] = mask[1] && mask[5];
1564 mask[8] = mask[0] && mask[4];
1566 mask[9] = mask[1] && mask[4];
1568 mask[10] = mask[0] && mask[2];
1570 mask[11] = mask[1] && mask[2];
1572 mask[12] = mask[0] && mask[3];
1574 mask[13] = mask[1] && mask[3];
1576 mask[14] = mask[3] && mask[4];
1578 mask[15] = mask[3] && mask[5];
1580 mask[16] = mask[2] && mask[4];
1582 mask[17] = mask[2] && mask[5];
1586 mask[18] = mask[1] && mask[3] && mask[5];
1588 mask[19] = mask[1] && mask[3] && mask[4];
1590 mask[20] = mask[0] && mask[3] && mask[4];
1592 mask[21] = mask[0] && mask[3] && mask[5];
1594 mask[22] = mask[1] && mask[2] && mask[5];
1596 mask[23] = mask[1] && mask[2] && mask[4];
1598 mask[24] = mask[0] && mask[2] && mask[4];
1600 mask[25] = mask[0] && mask[2] && mask[5];
1604 template<
typename Compare,
typename LeafNodeType>
1606 checkNeighbours(
const Index pos,
const typename LeafNodeType::ValueType * data,
bool (&mask)[26])
1608 using NodeT = LeafNodeType;
1611 if (mask[5] && Compare::check(data[pos - 1]))
return true;
1613 if (mask[4] && Compare::check(data[pos + 1]))
return true;
1615 if (mask[3] && Compare::check(data[pos - NodeT::DIM]))
return true;
1617 if (mask[2] && Compare::check(data[pos + NodeT::DIM]))
return true;
1619 if (mask[1] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM]))
return true;
1621 if (mask[0] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM]))
return true;
1623 if (mask[6] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM]))
return true;
1625 if (mask[7] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - 1]))
return true;
1627 if (mask[8] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + 1]))
return true;
1629 if (mask[9] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + 1]))
return true;
1631 if (mask[10] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM]))
return true;
1633 if (mask[11] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM]))
return true;
1635 if (mask[12] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM]))
return true;
1637 if (mask[13] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM]))
return true;
1639 if (mask[14] && Compare::check(data[pos - NodeT::DIM + 1]))
return true;
1641 if (mask[15] && Compare::check(data[pos - NodeT::DIM - 1]))
return true;
1643 if (mask[16] && Compare::check(data[pos + NodeT::DIM + 1]))
return true;
1645 if (mask[17] && Compare::check(data[pos + NodeT::DIM - 1]))
return true;
1647 if (mask[18] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM - 1]))
return true;
1649 if (mask[19] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM + 1]))
return true;
1651 if (mask[20] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM + 1]))
return true;
1653 if (mask[21] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM - 1]))
return true;
1655 if (mask[22] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM - 1]))
return true;
1657 if (mask[23] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM + 1]))
return true;
1659 if (mask[24] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM + 1]))
return true;
1661 if (mask[25] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM - 1]))
return true;
1667 template<
typename Compare,
typename AccessorType>
1669 checkNeighbours(
const Coord& ijk, AccessorType& acc,
bool (&mask)[26])
1671 for (
Int32 m = 0; m < 26; ++m) {
1681 template<
typename TreeType>
1682 struct ValidateIntersectingVoxels
1684 using ValueType =
typename TreeType::ValueType;
1685 using LeafNodeType =
typename TreeType::LeafNodeType;
1687 struct IsNegative {
static bool check(
const ValueType v) {
return v < ValueType(0.0); } };
1689 ValidateIntersectingVoxels(TreeType& tree, std::vector<LeafNodeType*>& nodes)
1691 , mNodes(nodes.empty() ? nullptr : &nodes[0])
1695 void operator()(
const tbb::blocked_range<size_t>& range)
const
1697 tree::ValueAccessor<const TreeType> acc(*mTree);
1698 bool neighbourMask[26];
1700 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1702 LeafNodeType& node = *mNodes[n];
1703 ValueType* data = node.buffer().data();
1705 typename LeafNodeType::ValueOnCIter it;
1706 for (it = node.cbeginValueOn(); it; ++it) {
1708 const Index pos = it.pos();
1710 ValueType& dist = data[pos];
1711 if (dist < 0.0 || dist > 0.75)
continue;
1714 maskNodeInternalNeighbours<LeafNodeType>(pos, neighbourMask);
1716 const bool hasNegativeNeighbour =
1717 checkNeighbours<IsNegative, LeafNodeType>(pos, data, neighbourMask) ||
1718 checkNeighbours<IsNegative>(node.offsetToGlobalCoord(pos), acc, neighbourMask);
1720 if (!hasNegativeNeighbour) {
1722 dist = ValueType(0.75) + Tolerance<ValueType>::epsilon();
1728 TreeType *
const mTree;
1729 LeafNodeType **
const mNodes;
1733 template<
typename TreeType>
1734 struct RemoveSelfIntersectingSurface
1736 using ValueType =
typename TreeType::ValueType;
1737 using LeafNodeType =
typename TreeType::LeafNodeType;
1740 struct Comp {
static bool check(
const ValueType v) {
return !(v > ValueType(0.75)); } };
1742 RemoveSelfIntersectingSurface(std::vector<LeafNodeType*>& nodes,
1743 TreeType& distTree, Int32TreeType& indexTree)
1744 : mNodes(nodes.empty() ? nullptr : &nodes[0])
1745 , mDistTree(&distTree)
1746 , mIndexTree(&indexTree)
1750 void operator()(
const tbb::blocked_range<size_t>& range)
const
1752 tree::ValueAccessor<const TreeType> distAcc(*mDistTree);
1753 tree::ValueAccessor<Int32TreeType> idxAcc(*mIndexTree);
1754 bool neighbourMask[26];
1756 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1758 LeafNodeType& distNode = *mNodes[n];
1759 ValueType* data = distNode.buffer().data();
1761 typename Int32TreeType::LeafNodeType* idxNode =
1762 idxAcc.probeLeaf(distNode.origin());
1764 typename LeafNodeType::ValueOnCIter it;
1765 for (it = distNode.cbeginValueOn(); it; ++it) {
1767 const Index pos = it.pos();
1769 if (!(data[pos] > 0.75))
continue;
1772 maskNodeInternalNeighbours<LeafNodeType>(pos, neighbourMask);
1774 const bool hasBoundaryNeighbour =
1775 checkNeighbours<Comp, LeafNodeType>(pos, data, neighbourMask) ||
1776 checkNeighbours<Comp>(distNode.offsetToGlobalCoord(pos),distAcc,neighbourMask);
1778 if (!hasBoundaryNeighbour) {
1779 distNode.setValueOff(pos);
1780 idxNode->setValueOff(pos);
1786 LeafNodeType * *
const mNodes;
1787 TreeType *
const mDistTree;
1788 Int32TreeType *
const mIndexTree;
1795 template<
typename NodeType>
1796 struct ReleaseChildNodes
1798 ReleaseChildNodes(NodeType ** nodes) : mNodes(nodes) {}
1800 void operator()(
const tbb::blocked_range<size_t>& range)
const {
1802 using NodeMaskType =
typename NodeType::NodeMaskType;
1804 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1805 const_cast<NodeMaskType&
>(mNodes[n]->getChildMask()).setOff();
1809 NodeType **
const mNodes;
1813 template<
typename TreeType>
1815 releaseLeafNodes(TreeType& tree)
1817 using RootNodeType =
typename TreeType::RootNodeType;
1818 using NodeChainType =
typename RootNodeType::NodeChainType;
1819 using InternalNodeType =
typename NodeChainType::template Get<1>;
1821 std::vector<InternalNodeType*> nodes;
1822 tree.getNodes(nodes);
1824 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
1825 ReleaseChildNodes<InternalNodeType>(nodes.empty() ?
nullptr : &nodes[0]));
1829 template<
typename TreeType>
1830 struct StealUniqueLeafNodes
1832 using LeafNodeType =
typename TreeType::LeafNodeType;
1834 StealUniqueLeafNodes(TreeType& lhsTree, TreeType& rhsTree,
1835 std::vector<LeafNodeType*>& overlappingNodes)
1836 : mLhsTree(&lhsTree)
1837 , mRhsTree(&rhsTree)
1838 , mNodes(&overlappingNodes)
1842 void operator()()
const {
1844 std::vector<LeafNodeType*> rhsLeafNodes;
1846 rhsLeafNodes.reserve(mRhsTree->leafCount());
1849 mRhsTree->stealNodes(rhsLeafNodes);
1851 tree::ValueAccessor<TreeType> acc(*mLhsTree);
1853 for (
size_t n = 0, N = rhsLeafNodes.size(); n < N; ++n) {
1854 if (!acc.probeLeaf(rhsLeafNodes[n]->origin())) {
1855 acc.addLeaf(rhsLeafNodes[n]);
1857 mNodes->push_back(rhsLeafNodes[n]);
1863 TreeType *
const mLhsTree;
1864 TreeType *
const mRhsTree;
1865 std::vector<LeafNodeType*> *
const mNodes;
1869 template<
typename DistTreeType,
typename IndexTreeType>
1871 combineData(DistTreeType& lhsDist, IndexTreeType& lhsIdx,
1872 DistTreeType& rhsDist, IndexTreeType& rhsIdx)
1874 using DistLeafNodeType =
typename DistTreeType::LeafNodeType;
1875 using IndexLeafNodeType =
typename IndexTreeType::LeafNodeType;
1877 std::vector<DistLeafNodeType*> overlappingDistNodes;
1878 std::vector<IndexLeafNodeType*> overlappingIdxNodes;
1881 tbb::task_group tasks;
1882 tasks.run(StealUniqueLeafNodes<DistTreeType>(lhsDist, rhsDist, overlappingDistNodes));
1883 tasks.run(StealUniqueLeafNodes<IndexTreeType>(lhsIdx, rhsIdx, overlappingIdxNodes));
1887 if (!overlappingDistNodes.empty() && !overlappingIdxNodes.empty()) {
1888 tbb::parallel_for(tbb::blocked_range<size_t>(0, overlappingDistNodes.size()),
1889 CombineLeafNodes<DistTreeType>(lhsDist, lhsIdx,
1890 &overlappingDistNodes[0], &overlappingIdxNodes[0]));
1900 template<
typename TreeType>
1901 struct VoxelizationData {
1903 using Ptr = std::unique_ptr<VoxelizationData>;
1904 using ValueType =
typename TreeType::ValueType;
1907 using UCharTreeType =
typename TreeType::template ValueConverter<unsigned char>::Type;
1909 using FloatTreeAcc = tree::ValueAccessor<TreeType>;
1910 using Int32TreeAcc = tree::ValueAccessor<Int32TreeType>;
1911 using UCharTreeAcc = tree::ValueAccessor<UCharTreeType>;
1915 : distTree(
std::numeric_limits<ValueType>::
max())
1918 , indexAcc(indexTree)
1919 , primIdTree(MaxPrimId)
1920 , primIdAcc(primIdTree)
1926 FloatTreeAcc distAcc;
1928 Int32TreeType indexTree;
1929 Int32TreeAcc indexAcc;
1931 UCharTreeType primIdTree;
1932 UCharTreeAcc primIdAcc;
1934 unsigned char getNewPrimId() {
1950 if (mPrimCount == MaxPrimId || primIdTree.leafCount() > 1000) {
1952 primIdTree.root().clear();
1953 primIdTree.clearAllAccessors();
1954 assert(mPrimCount == 0);
1957 return mPrimCount++;
1962 enum { MaxPrimId = 100 };
1964 unsigned char mPrimCount;
1968 template<
typename TreeType,
typename MeshDataAdapter,
typename Interrupter = util::NullInterrupter>
1969 class VoxelizePolygons
1973 using VoxelizationDataType = VoxelizationData<TreeType>;
1974 using DataTable = tbb::enumerable_thread_specific<typename VoxelizationDataType::Ptr>;
1976 VoxelizePolygons(DataTable& dataTable,
1978 Interrupter* interrupter =
nullptr)
1979 : mDataTable(&dataTable)
1981 , mInterrupter(interrupter)
1985 void operator()(
const tbb::blocked_range<size_t>& range)
const {
1987 typename VoxelizationDataType::Ptr& dataPtr = mDataTable->local();
1988 if (!dataPtr) dataPtr.reset(
new VoxelizationDataType());
1992 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
1995 thread::cancelGroupExecution();
1999 const size_t numVerts = mMesh->vertexCount(n);
2002 if (numVerts == 3 || numVerts == 4) {
2004 prim.index =
Int32(n);
2006 mMesh->getIndexSpacePoint(n, 0, prim.a);
2007 mMesh->getIndexSpacePoint(n, 1, prim.b);
2008 mMesh->getIndexSpacePoint(n, 2, prim.c);
2010 evalTriangle(prim, *dataPtr);
2012 if (numVerts == 4) {
2013 mMesh->getIndexSpacePoint(n, 3, prim.b);
2014 evalTriangle(prim, *dataPtr);
2022 bool wasInterrupted()
const {
return mInterrupter && mInterrupter->wasInterrupted(); }
2024 struct Triangle {
Vec3d a, b, c;
Int32 index; };
2028 enum { POLYGON_LIMIT = 1000 };
2030 SubTask(
const Triangle& prim, DataTable& dataTable,
2031 int subdivisionCount,
size_t polygonCount,
2032 Interrupter* interrupter =
nullptr)
2033 : mLocalDataTable(&dataTable)
2035 , mSubdivisionCount(subdivisionCount)
2036 , mPolygonCount(polygonCount)
2037 , mInterrupter(interrupter)
2041 void operator()()
const
2043 if (mSubdivisionCount <= 0 || mPolygonCount >= POLYGON_LIMIT) {
2045 typename VoxelizationDataType::Ptr& dataPtr = mLocalDataTable->local();
2046 if (!dataPtr) dataPtr.reset(
new VoxelizationDataType());
2048 voxelizeTriangle(mPrim, *dataPtr, mInterrupter);
2050 }
else if (!(mInterrupter && mInterrupter->wasInterrupted())) {
2051 spawnTasks(mPrim, *mLocalDataTable, mSubdivisionCount, mPolygonCount, mInterrupter);
2055 DataTable *
const mLocalDataTable;
2056 Triangle
const mPrim;
2057 int const mSubdivisionCount;
2058 size_t const mPolygonCount;
2059 Interrupter *
const mInterrupter;
2062 inline static int evalSubdivisionCount(
const Triangle& prim)
2064 const double ax = prim.a[0], bx = prim.b[0], cx = prim.c[0];
2067 const double ay = prim.a[1], by = prim.b[1], cy = prim.c[1];
2070 const double az = prim.a[2], bz = prim.b[2], cz = prim.c[2];
2073 return int(
std::max(dx,
std::max(dy, dz)) /
double(TreeType::LeafNodeType::DIM * 2));
2076 void evalTriangle(
const Triangle& prim, VoxelizationDataType& data)
const
2078 const size_t polygonCount = mMesh->polygonCount();
2079 const int subdivisionCount =
2080 polygonCount < SubTask::POLYGON_LIMIT ? evalSubdivisionCount(prim) : 0;
2082 if (subdivisionCount <= 0) {
2083 voxelizeTriangle(prim, data, mInterrupter);
2085 spawnTasks(prim, *mDataTable, subdivisionCount, polygonCount, mInterrupter);
2089 static void spawnTasks(
2090 const Triangle& mainPrim,
2091 DataTable& dataTable,
2092 int subdivisionCount,
2093 size_t polygonCount,
2094 Interrupter*
const interrupter)
2096 subdivisionCount -= 1;
2099 tbb::task_group tasks;
2101 const Vec3d ac = (mainPrim.a + mainPrim.c) * 0.5;
2102 const Vec3d bc = (mainPrim.b + mainPrim.c) * 0.5;
2103 const Vec3d ab = (mainPrim.a + mainPrim.b) * 0.5;
2106 prim.index = mainPrim.index;
2108 prim.a = mainPrim.a;
2111 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount, interrupter));
2116 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount, interrupter));
2119 prim.b = mainPrim.b;
2121 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount, interrupter));
2125 prim.c = mainPrim.c;
2126 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount, interrupter));
2131 static void voxelizeTriangle(
const Triangle& prim, VoxelizationDataType& data, Interrupter*
const interrupter)
2133 std::deque<Coord> coordList;
2137 coordList.push_back(ijk);
2142 updateDistance(ijk, prim, data);
2144 unsigned char primId = data.getNewPrimId();
2145 data.primIdAcc.setValueOnly(ijk, primId);
2147 while (!coordList.empty()) {
2148 if (interrupter && interrupter->wasInterrupted()) {
2149 thread::cancelGroupExecution();
2152 for (
Int32 pass = 0; pass < 1048576 && !coordList.empty(); ++pass) {
2153 ijk = coordList.back();
2154 coordList.pop_back();
2156 for (
Int32 i = 0; i < 26; ++i) {
2158 if (primId != data.primIdAcc.getValue(nijk)) {
2159 data.primIdAcc.setValueOnly(nijk, primId);
2160 if(updateDistance(nijk, prim, data)) coordList.push_back(nijk);
2167 static bool updateDistance(
const Coord& ijk,
const Triangle& prim, VoxelizationDataType& data)
2169 Vec3d uvw, voxelCenter(ijk[0], ijk[1], ijk[2]);
2171 using ValueType =
typename TreeType::ValueType;
2173 const ValueType dist = ValueType((voxelCenter -
2178 if (std::isnan(dist))
2181 const ValueType oldDist = data.distAcc.getValue(ijk);
2183 if (dist < oldDist) {
2184 data.distAcc.setValue(ijk, dist);
2185 data.indexAcc.setValue(ijk, prim.index);
2189 data.indexAcc.setValueOnly(ijk,
std::min(prim.index, data.indexAcc.getValue(ijk)));
2192 return !(dist > 0.75);
2195 DataTable *
const mDataTable;
2197 Interrupter *
const mInterrupter;
2204 template<
typename TreeType>
2205 struct DiffLeafNodeMask
2207 using AccessorType =
typename tree::ValueAccessor<TreeType>;
2208 using LeafNodeType =
typename TreeType::LeafNodeType;
2210 using BoolTreeType =
typename TreeType::template ValueConverter<bool>::Type;
2211 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
2213 DiffLeafNodeMask(
const TreeType& rhsTree,
2214 std::vector<BoolLeafNodeType*>& lhsNodes)
2215 : mRhsTree(&rhsTree), mLhsNodes(lhsNodes.empty() ? nullptr : &lhsNodes[0])
2219 void operator()(
const tbb::blocked_range<size_t>& range)
const {
2221 tree::ValueAccessor<const TreeType> acc(*mRhsTree);
2223 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2225 BoolLeafNodeType* lhsNode = mLhsNodes[n];
2226 const LeafNodeType* rhsNode = acc.probeConstLeaf(lhsNode->origin());
2228 if (rhsNode) lhsNode->topologyDifference(*rhsNode,
false);
2233 TreeType
const *
const mRhsTree;
2234 BoolLeafNodeType **
const mLhsNodes;
2238 template<
typename LeafNodeTypeA,
typename LeafNodeTypeB>
2239 struct UnionValueMasks
2241 UnionValueMasks(std::vector<LeafNodeTypeA*>& nodesA, std::vector<LeafNodeTypeB*>& nodesB)
2242 : mNodesA(nodesA.empty() ? nullptr : &nodesA[0])
2243 , mNodesB(nodesB.empty() ? nullptr : &nodesB[0])
2247 void operator()(
const tbb::blocked_range<size_t>& range)
const {
2248 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2249 mNodesA[n]->topologyUnion(*mNodesB[n]);
2254 LeafNodeTypeA **
const mNodesA;
2255 LeafNodeTypeB **
const mNodesB;
2259 template<
typename TreeType>
2260 struct ConstructVoxelMask
2262 using LeafNodeType =
typename TreeType::LeafNodeType;
2264 using BoolTreeType =
typename TreeType::template ValueConverter<bool>::Type;
2265 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
2267 ConstructVoxelMask(BoolTreeType& maskTree,
const TreeType& tree,
2268 std::vector<LeafNodeType*>& nodes)
2270 , mNodes(nodes.empty() ? nullptr : &nodes[0])
2271 , mLocalMaskTree(false)
2272 , mMaskTree(&maskTree)
2276 ConstructVoxelMask(ConstructVoxelMask& rhs, tbb::split)
2278 , mNodes(rhs.mNodes)
2279 , mLocalMaskTree(false)
2280 , mMaskTree(&mLocalMaskTree)
2284 void operator()(
const tbb::blocked_range<size_t>& range)
2286 using Iterator =
typename LeafNodeType::ValueOnCIter;
2288 tree::ValueAccessor<const TreeType> acc(*mTree);
2289 tree::ValueAccessor<BoolTreeType> maskAcc(*mMaskTree);
2291 Coord ijk, nijk, localCorod;
2294 for (
size_t n = range.begin(); n != range.end(); ++n) {
2296 LeafNodeType& node = *mNodes[n];
2298 CoordBBox bbox = node.getNodeBoundingBox();
2301 BoolLeafNodeType& maskNode = *maskAcc.touchLeaf(node.origin());
2303 for (Iterator it = node.cbeginValueOn(); it; ++it) {
2304 ijk = it.getCoord();
2307 localCorod = LeafNodeType::offsetToLocalCoord(pos);
2309 if (localCorod[2] <
int(LeafNodeType::DIM - 1)) {
2311 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2313 nijk = ijk.offsetBy(0, 0, 1);
2314 if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2317 if (localCorod[2] > 0) {
2319 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2321 nijk = ijk.offsetBy(0, 0, -1);
2322 if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2325 if (localCorod[1] <
int(LeafNodeType::DIM - 1)) {
2326 npos = pos + LeafNodeType::DIM;
2327 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2329 nijk = ijk.offsetBy(0, 1, 0);
2330 if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2333 if (localCorod[1] > 0) {
2334 npos = pos - LeafNodeType::DIM;
2335 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2337 nijk = ijk.offsetBy(0, -1, 0);
2338 if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2341 if (localCorod[0] <
int(LeafNodeType::DIM - 1)) {
2342 npos = pos + LeafNodeType::DIM * LeafNodeType::DIM;
2343 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2345 nijk = ijk.offsetBy(1, 0, 0);
2346 if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2349 if (localCorod[0] > 0) {
2350 npos = pos - LeafNodeType::DIM * LeafNodeType::DIM;
2351 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2353 nijk = ijk.offsetBy(-1, 0, 0);
2354 if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2360 void join(ConstructVoxelMask& rhs) { mMaskTree->merge(*rhs.mMaskTree); }
2363 TreeType
const *
const mTree;
2364 LeafNodeType **
const mNodes;
2366 BoolTreeType mLocalMaskTree;
2367 BoolTreeType *
const mMaskTree;
2372 template<
typename TreeType,
typename MeshDataAdapter>
2373 struct ExpandNarrowband
2375 using ValueType =
typename TreeType::ValueType;
2376 using LeafNodeType =
typename TreeType::LeafNodeType;
2377 using NodeMaskType =
typename LeafNodeType::NodeMaskType;
2379 using Int32LeafNodeType =
typename Int32TreeType::LeafNodeType;
2380 using BoolTreeType =
typename TreeType::template ValueConverter<bool>::Type;
2381 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
2388 Fragment() : idx(0), x(0), y(0), z(0), dist(0.0) {}
2391 : idx(idx_), x(x_), y(y_), z(z_), dist(dist_)
2395 bool operator<(
const Fragment& rhs)
const {
return idx < rhs.idx; }
2401 std::vector<BoolLeafNodeType*>& maskNodes,
2402 BoolTreeType& maskTree,
2404 Int32TreeType& indexTree,
2406 ValueType exteriorBandWidth,
2407 ValueType interiorBandWidth,
2408 ValueType voxelSize)
2409 : mMaskNodes(maskNodes.empty() ? nullptr : &maskNodes[0])
2410 , mMaskTree(&maskTree)
2411 , mDistTree(&distTree)
2412 , mIndexTree(&indexTree)
2414 , mNewMaskTree(false)
2416 , mUpdatedDistNodes()
2418 , mUpdatedIndexNodes()
2419 , mExteriorBandWidth(exteriorBandWidth)
2420 , mInteriorBandWidth(interiorBandWidth)
2421 , mVoxelSize(voxelSize)
2425 ExpandNarrowband(
const ExpandNarrowband& rhs, tbb::split)
2426 : mMaskNodes(rhs.mMaskNodes)
2427 , mMaskTree(rhs.mMaskTree)
2428 , mDistTree(rhs.mDistTree)
2429 , mIndexTree(rhs.mIndexTree)
2431 , mNewMaskTree(false)
2433 , mUpdatedDistNodes()
2435 , mUpdatedIndexNodes()
2436 , mExteriorBandWidth(rhs.mExteriorBandWidth)
2437 , mInteriorBandWidth(rhs.mInteriorBandWidth)
2438 , mVoxelSize(rhs.mVoxelSize)
2442 void join(ExpandNarrowband& rhs)
2444 mDistNodes.insert(mDistNodes.end(), rhs.mDistNodes.begin(), rhs.mDistNodes.end());
2445 mIndexNodes.insert(mIndexNodes.end(), rhs.mIndexNodes.begin(), rhs.mIndexNodes.end());
2447 mUpdatedDistNodes.insert(mUpdatedDistNodes.end(),
2448 rhs.mUpdatedDistNodes.begin(), rhs.mUpdatedDistNodes.end());
2450 mUpdatedIndexNodes.insert(mUpdatedIndexNodes.end(),
2451 rhs.mUpdatedIndexNodes.begin(), rhs.mUpdatedIndexNodes.end());
2453 mNewMaskTree.merge(rhs.mNewMaskTree);
2456 void operator()(
const tbb::blocked_range<size_t>& range)
2458 tree::ValueAccessor<BoolTreeType> newMaskAcc(mNewMaskTree);
2459 tree::ValueAccessor<TreeType> distAcc(*mDistTree);
2460 tree::ValueAccessor<Int32TreeType> indexAcc(*mIndexTree);
2462 std::vector<Fragment> fragments;
2463 fragments.reserve(256);
2465 std::unique_ptr<LeafNodeType> newDistNodePt;
2466 std::unique_ptr<Int32LeafNodeType> newIndexNodePt;
2468 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2470 BoolLeafNodeType& maskNode = *mMaskNodes[n];
2471 if (maskNode.isEmpty())
continue;
2475 const Coord& origin = maskNode.origin();
2477 LeafNodeType * distNodePt = distAcc.probeLeaf(origin);
2478 Int32LeafNodeType * indexNodePt = indexAcc.probeLeaf(origin);
2480 assert(!distNodePt == !indexNodePt);
2482 bool usingNewNodes =
false;
2484 if (!distNodePt && !indexNodePt) {
2486 const ValueType backgroundDist = distAcc.getValue(origin);
2488 if (!newDistNodePt.get() && !newIndexNodePt.get()) {
2489 newDistNodePt.reset(
new LeafNodeType(origin, backgroundDist));
2490 newIndexNodePt.reset(
new Int32LeafNodeType(origin, indexAcc.getValue(origin)));
2493 if ((backgroundDist < ValueType(0.0)) !=
2494 (newDistNodePt->getValue(0) < ValueType(0.0))) {
2495 newDistNodePt->buffer().fill(backgroundDist);
2498 newDistNodePt->setOrigin(origin);
2499 newIndexNodePt->setOrigin(origin);
2502 distNodePt = newDistNodePt.get();
2503 indexNodePt = newIndexNodePt.get();
2505 usingNewNodes =
true;
2512 for (
typename BoolLeafNodeType::ValueOnIter it = maskNode.beginValueOn(); it; ++it) {
2513 bbox.expand(it.getCoord());
2518 gatherFragments(fragments, bbox, distAcc, indexAcc);
2523 bbox = maskNode.getNodeBoundingBox();
2525 bool updatedLeafNodes =
false;
2527 for (
typename BoolLeafNodeType::ValueOnIter it = maskNode.beginValueOn(); it; ++it) {
2529 const Coord ijk = it.getCoord();
2531 if (updateVoxel(ijk, 5, fragments, *distNodePt, *indexNodePt, &updatedLeafNodes)) {
2533 for (
Int32 i = 0; i < 6; ++i) {
2535 if (bbox.isInside(nijk)) {
2536 mask.setOn(BoolLeafNodeType::coordToOffset(nijk));
2538 newMaskAcc.setValueOn(nijk);
2542 for (
Int32 i = 6; i < 26; ++i) {
2544 if (bbox.isInside(nijk)) {
2545 mask.setOn(BoolLeafNodeType::coordToOffset(nijk));
2551 if (updatedLeafNodes) {
2554 mask -= indexNodePt->getValueMask();
2556 for (
typename NodeMaskType::OnIterator it = mask.beginOn(); it; ++it) {
2558 const Index pos = it.pos();
2559 const Coord ijk = maskNode.origin() + LeafNodeType::offsetToLocalCoord(pos);
2561 if (updateVoxel(ijk, 6, fragments, *distNodePt, *indexNodePt)) {
2562 for (
Int32 i = 0; i < 6; ++i) {
2569 if (usingNewNodes) {
2570 newDistNodePt->topologyUnion(*newIndexNodePt);
2571 mDistNodes.push_back(newDistNodePt.release());
2572 mIndexNodes.push_back(newIndexNodePt.release());
2574 mUpdatedDistNodes.push_back(distNodePt);
2575 mUpdatedIndexNodes.push_back(indexNodePt);
2583 BoolTreeType& newMaskTree() {
return mNewMaskTree; }
2585 std::vector<LeafNodeType*>& newDistNodes() {
return mDistNodes; }
2586 std::vector<LeafNodeType*>& updatedDistNodes() {
return mUpdatedDistNodes; }
2588 std::vector<Int32LeafNodeType*>& newIndexNodes() {
return mIndexNodes; }
2589 std::vector<Int32LeafNodeType*>& updatedIndexNodes() {
return mUpdatedIndexNodes; }
2595 gatherFragments(std::vector<Fragment>& fragments,
const CoordBBox& bbox,
2596 tree::ValueAccessor<TreeType>& distAcc, tree::ValueAccessor<Int32TreeType>& indexAcc)
2599 const Coord nodeMin = bbox.min() & ~(LeafNodeType::DIM - 1);
2600 const Coord nodeMax = bbox.max() & ~(LeafNodeType::DIM - 1);
2605 for (ijk[0] = nodeMin[0]; ijk[0] <= nodeMax[0]; ijk[0] += LeafNodeType::DIM) {
2606 for (ijk[1] = nodeMin[1]; ijk[1] <= nodeMax[1]; ijk[1] += LeafNodeType::DIM) {
2607 for (ijk[2] = nodeMin[2]; ijk[2] <= nodeMax[2]; ijk[2] += LeafNodeType::DIM) {
2608 if (LeafNodeType* distleaf = distAcc.probeLeaf(ijk)) {
2611 ijk.offsetBy(LeafNodeType::DIM - 1));
2612 gatherFragments(fragments, region, *distleaf, *indexAcc.probeLeaf(ijk));
2618 std::sort(fragments.begin(), fragments.end());
2622 gatherFragments(std::vector<Fragment>& fragments,
const CoordBBox& bbox,
2623 const LeafNodeType& distLeaf,
const Int32LeafNodeType& idxLeaf)
const
2625 const typename LeafNodeType::NodeMaskType& mask = distLeaf.getValueMask();
2626 const ValueType* distData = distLeaf.buffer().data();
2627 const Int32* idxData = idxLeaf.buffer().data();
2629 for (
int x = bbox.min()[0]; x <= bbox.max()[0]; ++x) {
2630 const Index xPos = (x & (LeafNodeType::DIM - 1u)) << (2 * LeafNodeType::LOG2DIM);
2631 for (
int y = bbox.min()[1]; y <= bbox.max()[1]; ++y) {
2632 const Index yPos = xPos + ((y & (LeafNodeType::DIM - 1u)) << LeafNodeType::LOG2DIM);
2633 for (
int z = bbox.min()[2]; z <= bbox.max()[2]; ++z) {
2634 const Index pos = yPos + (z & (LeafNodeType::DIM - 1u));
2635 if (mask.isOn(pos)) {
2636 fragments.push_back(Fragment(idxData[pos],x,y,z, std::abs(distData[pos])));
2646 computeDistance(
const Coord& ijk,
const Int32 manhattanLimit,
2647 const std::vector<Fragment>& fragments,
Int32& closestPrimIdx)
const
2649 Vec3d a, b, c, uvw, voxelCenter(ijk[0], ijk[1], ijk[2]);
2653 for (
size_t n = 0, N = fragments.size(); n < N; ++n) {
2655 const Fragment& fragment = fragments[n];
2656 if (lastIdx == fragment.idx)
continue;
2658 const Int32 dx = std::abs(fragment.x - ijk[0]);
2659 const Int32 dy = std::abs(fragment.y - ijk[1]);
2660 const Int32 dz = std::abs(fragment.z - ijk[2]);
2662 const Int32 manhattan = dx + dy + dz;
2663 if (manhattan > manhattanLimit)
continue;
2665 lastIdx = fragment.idx;
2667 const size_t polygon = size_t(lastIdx);
2669 mMesh->getIndexSpacePoint(polygon, 0, a);
2670 mMesh->getIndexSpacePoint(polygon, 1, b);
2671 mMesh->getIndexSpacePoint(polygon, 2, c);
2673 primDist = (voxelCenter -
2677 if (4 == mMesh->vertexCount(polygon)) {
2679 mMesh->getIndexSpacePoint(polygon, 3, b);
2682 a, b, c, voxelCenter, uvw)).lengthSqr();
2684 if (tmpDist < primDist) primDist = tmpDist;
2687 if (primDist < dist) {
2689 closestPrimIdx = lastIdx;
2693 return ValueType(std::sqrt(dist)) * mVoxelSize;
2699 updateVoxel(
const Coord& ijk,
const Int32 manhattanLimit,
2700 const std::vector<Fragment>& fragments,
2701 LeafNodeType& distLeaf, Int32LeafNodeType& idxLeaf,
bool* updatedLeafNodes =
nullptr)
2703 Int32 closestPrimIdx = 0;
2704 const ValueType distance = computeDistance(ijk, manhattanLimit, fragments, closestPrimIdx);
2706 const Index pos = LeafNodeType::coordToOffset(ijk);
2707 const bool inside = distLeaf.getValue(pos) < ValueType(0.0);
2709 bool activateNeighbourVoxels =
false;
2711 if (!inside && distance < mExteriorBandWidth) {
2712 if (updatedLeafNodes) *updatedLeafNodes =
true;
2713 activateNeighbourVoxels = (distance + mVoxelSize) < mExteriorBandWidth;
2714 distLeaf.setValueOnly(pos, distance);
2715 idxLeaf.setValueOn(pos, closestPrimIdx);
2716 }
else if (inside && distance < mInteriorBandWidth) {
2717 if (updatedLeafNodes) *updatedLeafNodes =
true;
2718 activateNeighbourVoxels = (distance + mVoxelSize) < mInteriorBandWidth;
2719 distLeaf.setValueOnly(pos, -distance);
2720 idxLeaf.setValueOn(pos, closestPrimIdx);
2723 return activateNeighbourVoxels;
2728 BoolLeafNodeType **
const mMaskNodes;
2729 BoolTreeType *
const mMaskTree;
2730 TreeType *
const mDistTree;
2731 Int32TreeType *
const mIndexTree;
2735 BoolTreeType mNewMaskTree;
2737 std::vector<LeafNodeType*> mDistNodes, mUpdatedDistNodes;
2738 std::vector<Int32LeafNodeType*> mIndexNodes, mUpdatedIndexNodes;
2740 const ValueType mExteriorBandWidth, mInteriorBandWidth, mVoxelSize;
2744 template<
typename TreeType>
2746 using LeafNodeType =
typename TreeType::LeafNodeType;
2748 AddNodes(TreeType& tree, std::vector<LeafNodeType*>& nodes)
2749 : mTree(&tree) , mNodes(&nodes)
2753 void operator()()
const {
2754 tree::ValueAccessor<TreeType> acc(*mTree);
2755 std::vector<LeafNodeType*>& nodes = *mNodes;
2756 for (
size_t n = 0, N = nodes.size(); n < N; ++n) {
2757 acc.addLeaf(nodes[n]);
2761 TreeType *
const mTree;
2762 std::vector<LeafNodeType*> *
const mNodes;
2766 template<
typename TreeType,
typename Int32TreeType,
typename BoolTreeType,
typename MeshDataAdapter>
2770 Int32TreeType& indexTree,
2771 BoolTreeType& maskTree,
2772 std::vector<typename BoolTreeType::LeafNodeType*>& maskNodes,
2774 typename TreeType::ValueType exteriorBandWidth,
2775 typename TreeType::ValueType interiorBandWidth,
2776 typename TreeType::ValueType voxelSize)
2778 ExpandNarrowband<TreeType, MeshDataAdapter> expandOp(maskNodes, maskTree,
2779 distTree, indexTree, mesh, exteriorBandWidth, interiorBandWidth, voxelSize);
2781 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, maskNodes.size()), expandOp);
2783 tbb::parallel_for(tbb::blocked_range<size_t>(0, expandOp.updatedIndexNodes().size()),
2784 UnionValueMasks<typename TreeType::LeafNodeType, typename Int32TreeType::LeafNodeType>(
2785 expandOp.updatedDistNodes(), expandOp.updatedIndexNodes()));
2787 tbb::task_group tasks;
2788 tasks.run(AddNodes<TreeType>(distTree, expandOp.newDistNodes()));
2789 tasks.run(AddNodes<Int32TreeType>(indexTree, expandOp.newIndexNodes()));
2793 maskTree.merge(expandOp.newMaskTree());
2801 template<
typename TreeType>
2802 struct TransformValues
2804 using LeafNodeType =
typename TreeType::LeafNodeType;
2805 using ValueType =
typename TreeType::ValueType;
2807 TransformValues(std::vector<LeafNodeType*>& nodes,
2808 ValueType voxelSize,
bool unsignedDist)
2810 , mVoxelSize(voxelSize)
2811 , mUnsigned(unsignedDist)
2815 void operator()(
const tbb::blocked_range<size_t>& range)
const {
2817 typename LeafNodeType::ValueOnIter iter;
2819 const bool udf = mUnsigned;
2820 const ValueType w[2] = { -mVoxelSize, mVoxelSize };
2822 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2824 for (iter = mNodes[n]->beginValueOn(); iter; ++iter) {
2825 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
2826 val = w[udf || (val < ValueType(0.0))] * std::sqrt(std::abs(val));
2832 LeafNodeType * *
const mNodes;
2833 const ValueType mVoxelSize;
2834 const bool mUnsigned;
2839 template<
typename TreeType>
2840 struct InactivateValues
2842 using LeafNodeType =
typename TreeType::LeafNodeType;
2843 using ValueType =
typename TreeType::ValueType;
2845 InactivateValues(std::vector<LeafNodeType*>& nodes,
2846 ValueType exBandWidth, ValueType inBandWidth)
2847 : mNodes(nodes.empty() ? nullptr : &nodes[0])
2848 , mExBandWidth(exBandWidth)
2849 , mInBandWidth(inBandWidth)
2853 void operator()(
const tbb::blocked_range<size_t>& range)
const {
2855 typename LeafNodeType::ValueOnIter iter;
2856 const ValueType exVal = mExBandWidth;
2857 const ValueType inVal = -mInBandWidth;
2859 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2861 for (iter = mNodes[n]->beginValueOn(); iter; ++iter) {
2863 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
2865 const bool inside = val < ValueType(0.0);
2867 if (inside && !(val > inVal)) {
2870 }
else if (!inside && !(val < exVal)) {
2879 LeafNodeType * *
const mNodes;
2880 const ValueType mExBandWidth, mInBandWidth;
2884 template<
typename TreeType>
2887 using LeafNodeType =
typename TreeType::LeafNodeType;
2888 using ValueType =
typename TreeType::ValueType;
2890 OffsetValues(std::vector<LeafNodeType*>& nodes, ValueType offset)
2891 : mNodes(nodes.empty() ? nullptr : &nodes[0]), mOffset(offset)
2895 void operator()(
const tbb::blocked_range<size_t>& range)
const {
2897 const ValueType offset = mOffset;
2899 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2901 typename LeafNodeType::ValueOnIter iter = mNodes[n]->beginValueOn();
2903 for (; iter; ++iter) {
2904 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
2911 LeafNodeType * *
const mNodes;
2912 const ValueType mOffset;
2916 template<
typename TreeType>
2919 using LeafNodeType =
typename TreeType::LeafNodeType;
2920 using ValueType =
typename TreeType::ValueType;
2922 Renormalize(
const TreeType& tree,
const std::vector<LeafNodeType*>& nodes,
2923 ValueType* buffer, ValueType voxelSize)
2925 , mNodes(nodes.empty() ? nullptr : &nodes[0])
2927 , mVoxelSize(voxelSize)
2931 void operator()(
const tbb::blocked_range<size_t>& range)
const
2933 using Vec3Type = math::Vec3<ValueType>;
2935 tree::ValueAccessor<const TreeType> acc(*mTree);
2940 const ValueType dx = mVoxelSize, invDx = ValueType(1.0) / mVoxelSize;
2942 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2944 ValueType* bufferData = &mBuffer[n * LeafNodeType::SIZE];
2946 typename LeafNodeType::ValueOnCIter iter = mNodes[n]->cbeginValueOn();
2947 for (; iter; ++iter) {
2949 const ValueType phi0 = *iter;
2951 ijk = iter.getCoord();
2953 up[0] = acc.getValue(ijk.offsetBy(1, 0, 0)) - phi0;
2954 up[1] = acc.getValue(ijk.offsetBy(0, 1, 0)) - phi0;
2955 up[2] = acc.getValue(ijk.offsetBy(0, 0, 1)) - phi0;
2957 down[0] = phi0 - acc.getValue(ijk.offsetBy(-1, 0, 0));
2958 down[1] = phi0 - acc.getValue(ijk.offsetBy(0, -1, 0));
2959 down[2] = phi0 - acc.getValue(ijk.offsetBy(0, 0, -1));
2963 const ValueType diff =
math::Sqrt(normSqGradPhi) * invDx - ValueType(1.0);
2966 bufferData[iter.pos()] = phi0 - dx * S * diff;
2972 TreeType
const *
const mTree;
2973 LeafNodeType
const *
const *
const mNodes;
2974 ValueType *
const mBuffer;
2976 const ValueType mVoxelSize;
2980 template<
typename TreeType>
2983 using LeafNodeType =
typename TreeType::LeafNodeType;
2984 using ValueType =
typename TreeType::ValueType;
2986 MinCombine(std::vector<LeafNodeType*>& nodes,
const ValueType* buffer)
2987 : mNodes(nodes.empty() ? nullptr : &nodes[0]), mBuffer(buffer)
2991 void operator()(
const tbb::blocked_range<size_t>& range)
const {
2993 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
2995 const ValueType* bufferData = &mBuffer[n * LeafNodeType::SIZE];
2997 typename LeafNodeType::ValueOnIter iter = mNodes[n]->beginValueOn();
2999 for (; iter; ++iter) {
3000 ValueType& val =
const_cast<ValueType&
>(iter.getValue());
3001 val =
std::min(val, bufferData[iter.pos()]);
3007 LeafNodeType * *
const mNodes;
3008 ValueType
const *
const mBuffer;
3022 template <
typename FloatTreeT>
3026 using ConnectivityTable = mesh_to_volume_internal::LeafNodeConnectivityTable<FloatTreeT>;
3031 ConnectivityTable nodeConnectivity(tree);
3033 std::vector<size_t> zStartNodes, yStartNodes, xStartNodes;
3038 for (
size_t n = 0; n < nodeConnectivity.size(); ++n) {
3039 if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevX()[n]) {
3040 xStartNodes.push_back(n);
3043 if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevY()[n]) {
3044 yStartNodes.push_back(n);
3047 if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevZ()[n]) {
3048 zStartNodes.push_back(n);
3052 using SweepingOp = mesh_to_volume_internal::SweepExteriorSign<FloatTreeT>;
3057 tbb::parallel_for(tbb::blocked_range<size_t>(0, zStartNodes.size()),
3060 tbb::parallel_for(tbb::blocked_range<size_t>(0, yStartNodes.size()),
3063 tbb::parallel_for(tbb::blocked_range<size_t>(0, xStartNodes.size()),
3066 const size_t numLeafNodes = nodeConnectivity.size();
3067 const size_t numVoxels = numLeafNodes * FloatTreeT::LeafNodeType::SIZE;
3069 std::unique_ptr<bool[]> changedNodeMaskA{
new bool[numLeafNodes]};
3070 std::unique_ptr<bool[]> changedNodeMaskB{
new bool[numLeafNodes]};
3071 std::unique_ptr<bool[]> changedVoxelMask{
new bool[numVoxels]};
3073 mesh_to_volume_internal::fillArray(changedNodeMaskA.get(),
true, numLeafNodes);
3074 mesh_to_volume_internal::fillArray(changedNodeMaskB.get(),
false, numLeafNodes);
3075 mesh_to_volume_internal::fillArray(changedVoxelMask.get(),
false, numVoxels);
3077 const tbb::blocked_range<size_t> nodeRange(0, numLeafNodes);
3079 bool nodesUpdated =
false;
3084 tbb::parallel_for(nodeRange, mesh_to_volume_internal::SeedFillExteriorSign<FloatTreeT>(
3085 nodeConnectivity.nodes(), changedNodeMaskA.get()));
3093 tbb::parallel_for(nodeRange, mesh_to_volume_internal::SeedPoints<FloatTreeT>(
3094 nodeConnectivity, changedNodeMaskA.get(), changedNodeMaskB.get(),
3095 changedVoxelMask.get()));
3099 changedNodeMaskA.swap(changedNodeMaskB);
3101 nodesUpdated =
false;
3102 for (
size_t n = 0; n < numLeafNodes; ++n) {
3103 nodesUpdated |= changedNodeMaskA[n];
3104 if (nodesUpdated)
break;
3110 tbb::parallel_for(nodeRange, mesh_to_volume_internal::SyncVoxelMask<FloatTreeT>(
3111 nodeConnectivity.nodes(), changedNodeMaskA.get(), changedVoxelMask.get()));
3113 }
while (nodesUpdated);
3121 template <
typename Gr
idType,
typename MeshDataAdapter,
typename Interrupter>
3122 inline typename GridType::Ptr
3124 Interrupter& interrupter,
3127 float exteriorBandWidth,
3128 float interiorBandWidth,
3132 using GridTypePtr =
typename GridType::Ptr;
3133 using TreeType =
typename GridType::TreeType;
3134 using LeafNodeType =
typename TreeType::LeafNodeType;
3135 using ValueType =
typename GridType::ValueType;
3138 using Int32TreeType =
typename Int32GridType::TreeType;
3140 using BoolTreeType =
typename TreeType::template ValueConverter<bool>::Type;
3147 distGrid->setTransform(transform.
copy());
3149 ValueType exteriorWidth = ValueType(exteriorBandWidth);
3150 ValueType interiorWidth = ValueType(interiorBandWidth);
3154 if (!std::isfinite(exteriorWidth) || std::isnan(interiorWidth)) {
3155 std::stringstream msg;
3156 msg <<
"Illegal narrow band width: exterior = " << exteriorWidth
3157 <<
", interior = " << interiorWidth;
3162 const ValueType voxelSize = ValueType(transform.
voxelSize()[0]);
3164 if (!std::isfinite(voxelSize) ||
math::isZero(voxelSize)) {
3165 std::stringstream msg;
3166 msg <<
"Illegal transform, voxel size = " << voxelSize;
3172 exteriorWidth *= voxelSize;
3176 interiorWidth *= voxelSize;
3184 Int32GridType* indexGrid =
nullptr;
3186 typename Int32GridType::Ptr temporaryIndexGrid;
3188 if (polygonIndexGrid) {
3189 indexGrid = polygonIndexGrid;
3192 indexGrid = temporaryIndexGrid.get();
3195 indexGrid->newTree();
3196 indexGrid->setTransform(transform.
copy());
3198 if (computeSignedDistanceField) {
3202 interiorWidth = ValueType(0.0);
3205 TreeType& distTree = distGrid->tree();
3206 Int32TreeType& indexTree = indexGrid->tree();
3214 using VoxelizationDataType = mesh_to_volume_internal::VoxelizationData<TreeType>;
3215 using DataTable = tbb::enumerable_thread_specific<typename VoxelizationDataType::Ptr>;
3219 mesh_to_volume_internal::VoxelizePolygons<TreeType, MeshDataAdapter, Interrupter>;
3221 const tbb::blocked_range<size_t> polygonRange(0, mesh.polygonCount());
3223 tbb::parallel_for(polygonRange, Voxelizer(data, mesh, &interrupter));
3225 for (
typename DataTable::iterator i = data.begin(); i != data.end(); ++i) {
3226 VoxelizationDataType& dataItem = **i;
3227 mesh_to_volume_internal::combineData(
3228 distTree, indexTree, dataItem.distTree, dataItem.indexTree);
3234 if (interrupter.wasInterrupted(30))
return distGrid;
3241 if (computeSignedDistanceField) {
3246 std::vector<LeafNodeType*> nodes;
3247 nodes.reserve(distTree.leafCount());
3248 distTree.getNodes(nodes);
3250 const tbb::blocked_range<size_t> nodeRange(0, nodes.size());
3253 mesh_to_volume_internal::ComputeIntersectingVoxelSign<TreeType, MeshDataAdapter>;
3255 tbb::parallel_for(nodeRange, SignOp(nodes, distTree, indexTree, mesh));
3257 if (interrupter.wasInterrupted(45))
return distGrid;
3260 if (removeIntersectingVoxels) {
3262 tbb::parallel_for(nodeRange,
3263 mesh_to_volume_internal::ValidateIntersectingVoxels<TreeType>(distTree, nodes));
3265 tbb::parallel_for(nodeRange,
3266 mesh_to_volume_internal::RemoveSelfIntersectingSurface<TreeType>(
3267 nodes, distTree, indexTree));
3274 if (interrupter.wasInterrupted(50))
return distGrid;
3276 if (distTree.activeVoxelCount() == 0) {
3278 distTree.root().setBackground(exteriorWidth,
false);
3284 std::vector<LeafNodeType*> nodes;
3285 nodes.reserve(distTree.leafCount());
3286 distTree.getNodes(nodes);
3288 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3289 mesh_to_volume_internal::TransformValues<TreeType>(
3290 nodes, voxelSize, !computeSignedDistanceField));
3294 if (computeSignedDistanceField) {
3295 distTree.root().setBackground(exteriorWidth,
false);
3301 if (interrupter.wasInterrupted(54))
return distGrid;
3308 const ValueType minBandWidth = voxelSize * ValueType(2.0);
3310 if (interiorWidth > minBandWidth || exteriorWidth > minBandWidth) {
3313 BoolTreeType maskTree(
false);
3316 std::vector<LeafNodeType*> nodes;
3317 nodes.reserve(distTree.leafCount());
3318 distTree.getNodes(nodes);
3320 mesh_to_volume_internal::ConstructVoxelMask<TreeType> op(maskTree, distTree, nodes);
3321 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), op);
3327 float progress = 54.0f, step = 0.0f;
3329 2.0 * std::ceil((
std::max(interiorWidth, exteriorWidth) - minBandWidth) / voxelSize);
3331 if (estimated <
double(maxIterations)) {
3332 maxIterations = unsigned(estimated);
3333 step = 40.0f / float(maxIterations);
3336 std::vector<typename BoolTreeType::LeafNodeType*> maskNodes;
3341 if (interrupter.wasInterrupted(
int(progress)))
return distGrid;
3343 const size_t maskNodeCount = maskTree.leafCount();
3344 if (maskNodeCount == 0)
break;
3347 maskNodes.reserve(maskNodeCount);
3348 maskTree.getNodes(maskNodes);
3350 const tbb::blocked_range<size_t> range(0, maskNodes.size());
3352 tbb::parallel_for(range,
3353 mesh_to_volume_internal::DiffLeafNodeMask<TreeType>(distTree, maskNodes));
3355 mesh_to_volume_internal::expandNarrowband(distTree, indexTree, maskTree, maskNodes,
3356 mesh, exteriorWidth, interiorWidth, voxelSize);
3358 if ((++count) >= maxIterations)
break;
3363 if (interrupter.wasInterrupted(94))
return distGrid;
3365 if (!polygonIndexGrid) indexGrid->clear();
3373 if (computeSignedDistanceField && renormalizeValues) {
3375 std::vector<LeafNodeType*> nodes;
3376 nodes.reserve(distTree.leafCount());
3377 distTree.getNodes(nodes);
3379 std::unique_ptr<ValueType[]> buffer{
new ValueType[LeafNodeType::SIZE * nodes.size()]};
3381 const ValueType offset = ValueType(0.8 * voxelSize);
3383 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3384 mesh_to_volume_internal::OffsetValues<TreeType>(nodes, -offset));
3386 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3387 mesh_to_volume_internal::Renormalize<TreeType>(
3388 distTree, nodes, buffer.get(), voxelSize));
3390 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3391 mesh_to_volume_internal::MinCombine<TreeType>(nodes, buffer.get()));
3393 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3394 mesh_to_volume_internal::OffsetValues<TreeType>(
3395 nodes, offset - mesh_to_volume_internal::Tolerance<ValueType>::epsilon()));
3398 if (interrupter.wasInterrupted(99))
return distGrid;
3405 if (trimNarrowBand &&
std::min(interiorWidth, exteriorWidth) < voxelSize * ValueType(4.0)) {
3407 std::vector<LeafNodeType*> nodes;
3408 nodes.reserve(distTree.leafCount());
3409 distTree.getNodes(nodes);
3411 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3412 mesh_to_volume_internal::InactivateValues<TreeType>(
3413 nodes, exteriorWidth, computeSignedDistanceField ? interiorWidth : exteriorWidth));
3416 distTree, exteriorWidth, computeSignedDistanceField ? -interiorWidth : -exteriorWidth);
3423 template <
typename Gr
idType,
typename MeshDataAdapter>
3424 inline typename GridType::Ptr
3428 float exteriorBandWidth,
3429 float interiorBandWidth,
3434 return meshToVolume<GridType>(nullInterrupter, mesh, transform,
3435 exteriorBandWidth, interiorBandWidth, flags, polygonIndexGrid);
3446 template<
typename Gr
idType,
typename Interrupter>
3447 inline typename std::enable_if<std::is_floating_point<typename GridType::ValueType>::value,
3448 typename GridType::Ptr>::type
3450 Interrupter& interrupter,
3451 const openvdb::math::Transform& xform,
3452 const std::vector<Vec3s>& points,
3453 const std::vector<Vec3I>& triangles,
3454 const std::vector<Vec4I>& quads,
3457 bool unsignedDistanceField =
false)
3459 if (points.empty()) {
3460 return typename GridType::Ptr(
new GridType(
typename GridType::ValueType(exBandWidth)));
3463 const size_t numPoints = points.size();
3464 std::unique_ptr<Vec3s[]> indexSpacePoints{
new Vec3s[numPoints]};
3467 tbb::parallel_for(tbb::blocked_range<size_t>(0, numPoints),
3468 mesh_to_volume_internal::TransformPoints<Vec3s>(
3469 &points[0], indexSpacePoints.get(), xform));
3473 if (quads.empty()) {
3475 QuadAndTriangleDataAdapter<Vec3s, Vec3I>
3476 mesh(indexSpacePoints.get(), numPoints, &triangles[0], triangles.size());
3478 return meshToVolume<GridType>(
3479 interrupter, mesh, xform, exBandWidth, inBandWidth, conversionFlags);
3481 }
else if (triangles.empty()) {
3483 QuadAndTriangleDataAdapter<Vec3s, Vec4I>
3484 mesh(indexSpacePoints.get(), numPoints, &quads[0], quads.size());
3486 return meshToVolume<GridType>(
3487 interrupter, mesh, xform, exBandWidth, inBandWidth, conversionFlags);
3492 const size_t numPrimitives = triangles.size() + quads.size();
3493 std::unique_ptr<Vec4I[]> prims{
new Vec4I[numPrimitives]};
3495 for (
size_t n = 0, N = triangles.size(); n < N; ++n) {
3496 const Vec3I& triangle = triangles[n];
3497 Vec4I& prim = prims[n];
3498 prim[0] = triangle[0];
3499 prim[1] = triangle[1];
3500 prim[2] = triangle[2];
3504 const size_t offset = triangles.size();
3505 for (
size_t n = 0, N = quads.size(); n < N; ++n) {
3506 prims[offset + n] = quads[n];
3509 QuadAndTriangleDataAdapter<Vec3s, Vec4I>
3510 mesh(indexSpacePoints.get(), numPoints, prims.get(), numPrimitives);
3512 return meshToVolume<GridType>(interrupter, mesh, xform,
3513 exBandWidth, inBandWidth, conversionFlags);
3519 template<
typename Gr
idType,
typename Interrupter>
3520 inline typename std::enable_if<!std::is_floating_point<typename GridType::ValueType>::value,
3521 typename GridType::Ptr>::type
3524 const math::Transform& ,
3525 const std::vector<Vec3s>& ,
3526 const std::vector<Vec3I>& ,
3527 const std::vector<Vec4I>& ,
3533 "mesh to volume conversion is supported only for scalar floating-point grids");
3543 template<
typename Gr
idType>
3544 inline typename GridType::Ptr
3546 const openvdb::math::Transform& xform,
3547 const std::vector<Vec3s>& points,
3548 const std::vector<Vec3I>& triangles,
3552 std::vector<Vec4I> quads(0);
3553 return doMeshConversion<GridType>(nullInterrupter, xform, points, triangles, quads,
3554 halfWidth, halfWidth);
3558 template<
typename Gr
idType,
typename Interrupter>
3559 inline typename GridType::Ptr
3561 Interrupter& interrupter,
3562 const openvdb::math::Transform& xform,
3563 const std::vector<Vec3s>& points,
3564 const std::vector<Vec3I>& triangles,
3567 std::vector<Vec4I> quads(0);
3568 return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3569 halfWidth, halfWidth);
3573 template<
typename Gr
idType>
3574 inline typename GridType::Ptr
3576 const openvdb::math::Transform& xform,
3577 const std::vector<Vec3s>& points,
3578 const std::vector<Vec4I>& quads,
3582 std::vector<Vec3I> triangles(0);
3583 return doMeshConversion<GridType>(nullInterrupter, xform, points, triangles, quads,
3584 halfWidth, halfWidth);
3588 template<
typename Gr
idType,
typename Interrupter>
3589 inline typename GridType::Ptr
3591 Interrupter& interrupter,
3592 const openvdb::math::Transform& xform,
3593 const std::vector<Vec3s>& points,
3594 const std::vector<Vec4I>& quads,
3597 std::vector<Vec3I> triangles(0);
3598 return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3599 halfWidth, halfWidth);
3603 template<
typename Gr
idType>
3604 inline typename GridType::Ptr
3606 const openvdb::math::Transform& xform,
3607 const std::vector<Vec3s>& points,
3608 const std::vector<Vec3I>& triangles,
3609 const std::vector<Vec4I>& quads,
3613 return doMeshConversion<GridType>(nullInterrupter, xform, points, triangles, quads,
3614 halfWidth, halfWidth);
3618 template<
typename Gr
idType,
typename Interrupter>
3619 inline typename GridType::Ptr
3621 Interrupter& interrupter,
3622 const openvdb::math::Transform& xform,
3623 const std::vector<Vec3s>& points,
3624 const std::vector<Vec3I>& triangles,
3625 const std::vector<Vec4I>& quads,
3628 return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3629 halfWidth, halfWidth);
3633 template<
typename Gr
idType>
3634 inline typename GridType::Ptr
3636 const openvdb::math::Transform& xform,
3637 const std::vector<Vec3s>& points,
3638 const std::vector<Vec3I>& triangles,
3639 const std::vector<Vec4I>& quads,
3644 return doMeshConversion<GridType>(nullInterrupter, xform, points, triangles,
3645 quads, exBandWidth, inBandWidth);
3649 template<
typename Gr
idType,
typename Interrupter>
3650 inline typename GridType::Ptr
3652 Interrupter& interrupter,
3653 const openvdb::math::Transform& xform,
3654 const std::vector<Vec3s>& points,
3655 const std::vector<Vec3I>& triangles,
3656 const std::vector<Vec4I>& quads,
3660 return doMeshConversion<GridType>(interrupter, xform, points, triangles,
3661 quads, exBandWidth, inBandWidth);
3665 template<
typename Gr
idType>
3666 inline typename GridType::Ptr
3668 const openvdb::math::Transform& xform,
3669 const std::vector<Vec3s>& points,
3670 const std::vector<Vec3I>& triangles,
3671 const std::vector<Vec4I>& quads,
3675 return doMeshConversion<GridType>(nullInterrupter, xform, points, triangles, quads,
3676 bandWidth, bandWidth,
true);
3680 template<
typename Gr
idType,
typename Interrupter>
3681 inline typename GridType::Ptr
3683 Interrupter& interrupter,
3684 const openvdb::math::Transform& xform,
3685 const std::vector<Vec3s>& points,
3686 const std::vector<Vec3I>& triangles,
3687 const std::vector<Vec4I>& quads,
3690 return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3691 bandWidth, bandWidth,
true);
3699 inline std::ostream&
3702 ostr <<
"{[ " << rhs.
mXPrim <<
", " << rhs.
mXDist <<
"]";
3703 ostr <<
" [ " << rhs.
mYPrim <<
", " << rhs.
mYDist <<
"]";
3704 ostr <<
" [ " << rhs.
mZPrim <<
", " << rhs.
mZDist <<
"]}";
3709 inline MeshToVoxelEdgeData::EdgeData
3724 const std::vector<Vec3s>& pointList,
3725 const std::vector<Vec4I>& polygonList);
3727 void run(
bool threaded =
true);
3730 inline void operator() (
const tbb::blocked_range<size_t> &range);
3738 struct Primitive {
Vec3d a, b, c, d;
Int32 index; };
3740 template<
bool IsQuad>
3741 inline void voxelize(
const Primitive&);
3743 template<
bool IsQuad>
3744 inline bool evalPrimitive(
const Coord&,
const Primitive&);
3746 inline bool rayTriangleIntersection(
const Vec3d& origin,
const Vec3d& dir,
3753 const std::vector<Vec3s>& mPointList;
3754 const std::vector<Vec4I>& mPolygonList;
3757 using IntTreeT = TreeType::ValueConverter<Int32>::Type;
3758 IntTreeT mLastPrimTree;
3764 MeshToVoxelEdgeData::GenEdgeData::GenEdgeData(
3765 const std::vector<Vec3s>& pointList,
3766 const std::vector<Vec4I>& polygonList)
3769 , mPointList(pointList)
3770 , mPolygonList(polygonList)
3772 , mLastPrimAccessor(mLastPrimTree)
3781 , mPointList(rhs.mPointList)
3782 , mPolygonList(rhs.mPolygonList)
3784 , mLastPrimAccessor(mLastPrimTree)
3793 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPolygonList.size()), *
this);
3795 (*this)(tbb::blocked_range<size_t>(0, mPolygonList.size()));
3804 using NodeChainType = RootNodeType::NodeChainType;
3805 static_assert(NodeChainType::Size > 1,
"expected tree height > 1");
3806 using InternalNodeType =
typename NodeChainType::template Get<1>;
3814 for ( ; leafIt; ++leafIt) {
3815 ijk = leafIt->origin();
3821 mAccessor.addLeaf(rhs.mAccessor.
probeLeaf(ijk));
3822 InternalNodeType* node = rhs.mAccessor.
getNode<InternalNodeType>();
3824 rhs.mAccessor.
clear();
3834 if (!lhsLeafPt->isValueOn(offset)) {
3835 lhsLeafPt->setValueOn(offset, rhsValue);
3867 for (
size_t n = range.begin(); n < range.end(); ++n) {
3869 const Vec4I& verts = mPolygonList[n];
3871 prim.index =
Int32(n);
3872 prim.a =
Vec3d(mPointList[verts[0]]);
3873 prim.b =
Vec3d(mPointList[verts[1]]);
3874 prim.c =
Vec3d(mPointList[verts[2]]);
3877 prim.d =
Vec3d(mPointList[verts[3]]);
3878 voxelize<true>(prim);
3880 voxelize<false>(prim);
3886 template<
bool IsQuad>
3888 MeshToVoxelEdgeData::GenEdgeData::voxelize(
const Primitive& prim)
3890 std::deque<Coord> coordList;
3894 coordList.push_back(ijk);
3896 evalPrimitive<IsQuad>(ijk, prim);
3898 while (!coordList.empty()) {
3900 ijk = coordList.back();
3901 coordList.pop_back();
3903 for (
Int32 i = 0; i < 26; ++i) {
3906 if (prim.index != mLastPrimAccessor.getValue(nijk)) {
3907 mLastPrimAccessor.setValue(nijk, prim.index);
3908 if(evalPrimitive<IsQuad>(nijk, prim)) coordList.push_back(nijk);
3915 template<
bool IsQuad>
3917 MeshToVoxelEdgeData::GenEdgeData::evalPrimitive(
const Coord& ijk,
const Primitive& prim)
3919 Vec3d uvw, org(ijk[0], ijk[1], ijk[2]);
3920 bool intersecting =
false;
3924 mAccessor.probeValue(ijk, edgeData);
3927 double dist = (org -
3930 if (rayTriangleIntersection(org,
Vec3d(1.0, 0.0, 0.0), prim.a, prim.c, prim.b, t)) {
3931 if (t < edgeData.mXDist) {
3932 edgeData.mXDist = float(t);
3933 edgeData.mXPrim = prim.index;
3934 intersecting =
true;
3938 if (rayTriangleIntersection(org,
Vec3d(0.0, 1.0, 0.0), prim.a, prim.c, prim.b, t)) {
3939 if (t < edgeData.mYDist) {
3940 edgeData.mYDist = float(t);
3941 edgeData.mYPrim = prim.index;
3942 intersecting =
true;
3946 if (rayTriangleIntersection(org,
Vec3d(0.0, 0.0, 1.0), prim.a, prim.c, prim.b, t)) {
3947 if (t < edgeData.mZDist) {
3948 edgeData.mZDist = float(t);
3949 edgeData.mZPrim = prim.index;
3950 intersecting =
true;
3956 double secondDist = (org -
3959 if (secondDist < dist) dist = secondDist;
3961 if (rayTriangleIntersection(org,
Vec3d(1.0, 0.0, 0.0), prim.a, prim.d, prim.c, t)) {
3962 if (t < edgeData.mXDist) {
3963 edgeData.mXDist = float(t);
3964 edgeData.mXPrim = prim.index;
3965 intersecting =
true;
3969 if (rayTriangleIntersection(org,
Vec3d(0.0, 1.0, 0.0), prim.a, prim.d, prim.c, t)) {
3970 if (t < edgeData.mYDist) {
3971 edgeData.mYDist = float(t);
3972 edgeData.mYPrim = prim.index;
3973 intersecting =
true;
3977 if (rayTriangleIntersection(org,
Vec3d(0.0, 0.0, 1.0), prim.a, prim.d, prim.c, t)) {
3978 if (t < edgeData.mZDist) {
3979 edgeData.mZDist = float(t);
3980 edgeData.mZPrim = prim.index;
3981 intersecting =
true;
3986 if (intersecting) mAccessor.setValue(ijk, edgeData);
3988 return (dist < 0.86602540378443861);
3993 MeshToVoxelEdgeData::GenEdgeData::rayTriangleIntersection(
4004 double divisor = s1.
dot(e1);
4005 if (!(std::abs(divisor) > 0.0))
return false;
4009 double inv_divisor = 1.0 / divisor;
4010 Vec3d d = origin - a;
4011 double b1 = d.
dot(s1) * inv_divisor;
4013 if (b1 < 0.0 || b1 > 1.0)
return false;
4016 double b2 = dir.
dot(s2) * inv_divisor;
4018 if (b2 < 0.0 || (b1 + b2) > 1.0)
return false;
4022 t = e2.dot(s2) * inv_divisor;
4023 return (t < 0.0) ? false :
true;
4039 const std::vector<Vec3s>& pointList,
4040 const std::vector<Vec4I>& polygonList)
4054 std::vector<Vec3d>& points,
4055 std::vector<Index32>& primitives)
4065 point[0] = double(coord[0]) + data.
mXDist;
4066 point[1] = double(coord[1]);
4067 point[2] = double(coord[2]);
4069 points.push_back(point);
4070 primitives.push_back(data.
mXPrim);
4074 point[0] = double(coord[0]);
4075 point[1] = double(coord[1]) + data.
mYDist;
4076 point[2] = double(coord[2]);
4078 points.push_back(point);
4079 primitives.push_back(data.
mYPrim);
4083 point[0] = double(coord[0]);
4084 point[1] = double(coord[1]);
4085 point[2] = double(coord[2]) + data.
mZDist;
4087 points.push_back(point);
4088 primitives.push_back(data.
mZPrim);
4098 point[0] = double(coord[0]);
4099 point[1] = double(coord[1]) + data.
mYDist;
4100 point[2] = double(coord[2]);
4102 points.push_back(point);
4103 primitives.push_back(data.
mYPrim);
4107 point[0] = double(coord[0]);
4108 point[1] = double(coord[1]);
4109 point[2] = double(coord[2]) + data.
mZDist;
4111 points.push_back(point);
4112 primitives.push_back(data.
mZPrim);
4120 point[0] = double(coord[0]);
4121 point[1] = double(coord[1]) + data.
mYDist;
4122 point[2] = double(coord[2]);
4124 points.push_back(point);
4125 primitives.push_back(data.
mYPrim);
4134 point[0] = double(coord[0]) + data.
mXDist;
4135 point[1] = double(coord[1]);
4136 point[2] = double(coord[2]);
4138 points.push_back(point);
4139 primitives.push_back(data.
mXPrim);
4143 point[0] = double(coord[0]);
4144 point[1] = double(coord[1]) + data.
mYDist;
4145 point[2] = double(coord[2]);
4147 points.push_back(point);
4148 primitives.push_back(data.
mYPrim);
4158 point[0] = double(coord[0]) + data.
mXDist;
4159 point[1] = double(coord[1]);
4160 point[2] = double(coord[2]);
4162 points.push_back(point);
4163 primitives.push_back(data.
mXPrim);
4172 point[0] = double(coord[0]) + data.
mXDist;
4173 point[1] = double(coord[1]);
4174 point[2] = double(coord[2]);
4176 points.push_back(point);
4177 primitives.push_back(data.
mXPrim);
4181 point[0] = double(coord[0]);
4182 point[1] = double(coord[1]);
4183 point[2] = double(coord[2]) + data.
mZDist;
4185 points.push_back(point);
4186 primitives.push_back(data.
mZPrim);
4195 point[0] = double(coord[0]);
4196 point[1] = double(coord[1]);
4197 point[2] = double(coord[2]) + data.
mZDist;
4199 points.push_back(point);
4200 primitives.push_back(data.
mZPrim);
4206 template<
typename Gr
idType,
typename VecType>
4207 inline typename GridType::Ptr
4209 const openvdb::math::Transform& xform,
4210 typename VecType::ValueType halfWidth)
4216 points[0] =
Vec3s(pmin[0], pmin[1], pmin[2]);
4217 points[1] =
Vec3s(pmin[0], pmin[1], pmax[2]);
4218 points[2] =
Vec3s(pmax[0], pmin[1], pmax[2]);
4219 points[3] =
Vec3s(pmax[0], pmin[1], pmin[2]);
4220 points[4] =
Vec3s(pmin[0], pmax[1], pmin[2]);
4221 points[5] =
Vec3s(pmin[0], pmax[1], pmax[2]);
4222 points[6] =
Vec3s(pmax[0], pmax[1], pmax[2]);
4223 points[7] =
Vec3s(pmax[0], pmax[1], pmin[2]);
4226 faces[0] =
Vec4I(0, 1, 2, 3);
4227 faces[1] =
Vec4I(7, 6, 5, 4);
4228 faces[2] =
Vec4I(4, 5, 1, 0);
4229 faces[3] =
Vec4I(6, 7, 3, 2);
4230 faces[4] =
Vec4I(0, 3, 7, 4);
4231 faces[5] =
Vec4I(1, 5, 6, 2);
4235 return meshToVolume<GridType>(mesh, xform, halfWidth, halfWidth);
Efficient multi-threaded replacement of the background values in tree.
Defined various multi-threaded utility functions for trees.
Propagate the signs of distance values from the active voxels in the narrow band to the inactive valu...
Axis-aligned bounding box.
Definition: BBox.h:24
const Vec3T & max() const
Return a const reference to the maximum point of this bounding box.
Definition: BBox.h:64
const Vec3T & min() const
Return a const reference to the minimum point of this bounding box.
Definition: BBox.h:62
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:25
void minComponent(const Coord &other)
Perform a component-wise minimum with the other Coord.
Definition: Coord.h:175
void maxComponent(const Coord &other)
Perform a component-wise maximum with the other Coord.
Definition: Coord.h:183
static Coord min()
Return the smallest possible coordinate.
Definition: Coord.h:43
static Coord floor(const Vec3< T > &xyz)
Return the largest integer coordinates that are not greater than xyz (node centered conversion).
Definition: Coord.h:56
static Coord max()
Return the largest possible coordinate.
Definition: Coord.h:46
T dot(const Vec3< T > &v) const
Dot product.
Definition: Vec3.h:195
Vec3< T > cross(const Vec3< T > &v) const
Return the cross product of "this" vector and v;.
Definition: Vec3.h:224
Base class for tree-traversal iterators over all leaf nodes (but not leaf voxels)
Definition: TreeIterator.h:1187
Base class for tree-traversal iterators over tile and voxel values.
Definition: TreeIterator.h:617
const ValueT & getValue() const
Return the tile or voxel value to which this iterator is currently pointing.
Definition: TreeIterator.h:692
void clearAllAccessors()
Clear all registered accessors.
Definition: Tree.h:1377
void merge(Tree &other, MergePolicy=MERGE_ACTIVE_STATES)
Efficiently merge another tree into this tree using one of several schemes.
Definition: Tree.h:1690
_RootNodeType RootNodeType
Definition: Tree.h:181
LeafNodeType * probeLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains voxel (x, y, z). If no such node exists,...
Definition: Tree.h:1564
void clear()
Remove all tiles from this tree and all nodes other than the root node.
Definition: Tree.h:1318
typename RootNodeType::LeafNodeType LeafNodeType
Definition: Tree.h:184
LeafIter beginLeaf()
Return an iterator over all leaf nodes in this tree.
Definition: Tree.h:1015
LeafNodeT * probeLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains voxel (x, y, z), or nullptr if no such node exists.
Definition: ValueAccessor.h:379
bool probeValue(const Coord &xyz, ValueType &value) const
Return the active state of the voxel as well as its value.
Definition: ValueAccessor.h:229
NodeType * getNode()
Return the cached node of type NodeType. [Mainly for internal use].
Definition: ValueAccessor.h:304
void clear() override
Remove all nodes from this cache, then reinsert the root node.
Definition: ValueAccessor.h:393
Convert polygonal meshes that consist of quads and/or triangles into signed or unsigned distance fiel...
Partitions points into BucketLog2Dim aligned buckets using a parallel radix-based sorting algorithm.
#define OPENVDB_LOG_DEBUG(message)
In debug builds only, log a debugging message of the form 'someVar << "text" << .....
Definition: logging.h:266
void run(const char *ax, openvdb::GridBase &grid)
Run a full AX pipeline (parse, compile and execute) on a single OpenVDB Grid.
bool operator<(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition: Tuple.h:189
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:338
OPENVDB_API Vec3d closestPointOnTriangleToPoint(const Vec3d &a, const Vec3d &b, const Vec3d &c, const Vec3d &p, Vec3d &uvw)
Closest Point on Triangle to Point. Given a triangle abc and a point p, return the point on abc close...
bool operator>(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition: Tuple.h:201
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:764
Vec3< double > Vec3d
Definition: Vec3.h:668
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:326
Type Pow2(Type x)
Return x2.
Definition: Math.h:551
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:444
Axis
Definition: Math.h:904
@ Z_AXIS
Definition: Math.h:907
@ X_AXIS
Definition: Math.h:905
@ Y_AXIS
Definition: Math.h:906
Vec3< float > Vec3s
Definition: Vec3.h:667
OPENVDB_API const Coord COORD_OFFSETS[26]
coordinate offset table for neighboring voxels
OPENVDB_API const Index32 INVALID_IDX
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:49
static const Real LEVEL_SET_HALF_WIDTH
Definition: Types.h:343
Index32 Index
Definition: Types.h:54
math::Vec4< Index32 > Vec4I
Definition: Types.h:88
@ GRID_LEVEL_SET
Definition: Types.h:337
@ GRID_UNKNOWN
Definition: Types.h:336
uint32_t Index32
Definition: Types.h:52
math::Vec3< Index32 > Vec3I
Definition: Types.h:73
int32_t Int32
Definition: Types.h:56
Definition: Exceptions.h:13
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
Tree< typename RootNodeType::template ValueConverter< Int32 >::Type > Type
Definition: Tree.h:196
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:26
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:116
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:180