16 #ifndef OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
17 #define OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
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
115 float exteriorBandWidth = 3.0
f,
116 float interiorBandWidth = 3.0
f,
136 template <
typename Gr
idType,
typename MeshDataAdapter,
typename Interrupter>
137 inline typename GridType::Ptr
139 Interrupter& interrupter,
142 float exteriorBandWidth = 3.0
f,
143 float interiorBandWidth = 3.0
f,
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)
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,
237 float halfWidth =
float(LEVEL_SET_HALF_WIDTH));
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,
247 float halfWidth =
float(LEVEL_SET_HALF_WIDTH));
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,
271 float halfWidth =
float(LEVEL_SET_HALF_WIDTH));
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,
281 float halfWidth =
float(LEVEL_SET_HALF_WIDTH));
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,
307 float halfWidth =
float(LEVEL_SET_HALF_WIDTH));
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,
318 float halfWidth =
float(LEVEL_SET_HALF_WIDTH));
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,
410 typename VecType::ValueType halfWidth = LEVEL_SET_HALF_WIDTH);
422 template <
typename FloatTreeT>
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,
516 const math::Transform& xform)
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];
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();
603 const DistValueType& lhsValue = lhsDistData[
offset];
604 const DistValueType& rhsValue = rhsDistData[
offset];
606 if (rhsValue < lhsValue) {
607 lhsDistNode->setValueOn(
offset, rhsValue);
610 lhsIdxNode->setValueOn(
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()]};
773 StashOriginAndStoreOffset<TreeType>(mLeafNodes, coordinates.get()));
776 mOffsets.reset(
new size_t[mLeafNodes.size() * 6]);
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();
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();
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;
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]) {
1165 ValueType* data = mNodes[
n]->buffer().data();
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]) {
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]) {
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]) {
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);
1450 direction = nxyz - point;
1451 direction.normalize();
1454 dir1 = xyz - points[pointIndex];
1457 if (points[pointIndex + 1].
dot(dir1) > 0.0) {
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) {
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;
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])
1672 if (!mask[
m] &&
Compare::check(acc.getValue(ijk + util::COORD_OFFSETS[
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);
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()) {
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;
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) {
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(); }
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];
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()) {
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) {
2157 nijk = ijk + util::COORD_OFFSETS[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;
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;
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;
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) {
2534 const Coord nijk = ijk + util::COORD_OFFSETS[i];
2535 if (bbox.isInside(nijk)) {
2536 mask.setOn(BoolLeafNodeType::coordToOffset(nijk));
2538 newMaskAcc.setValueOn(nijk);
2542 for (
Int32 i = 6; i < 26; ++i) {
2543 const Coord nijk = ijk + util::COORD_OFFSETS[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) {
2563 newMaskAcc.setValueOn(ijk + util::COORD_OFFSETS[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;
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());
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) {
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) {
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>;
3066 const size_t numLeafNodes = nodeConnectivity.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()));
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;
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());
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>;
3257 if (interrupter.wasInterrupted(45))
return distGrid;
3260 if (removeIntersectingVoxels) {
3263 mesh_to_volume_internal::ValidateIntersectingVoxels<TreeType>(distTree, nodes));
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);
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());
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);
3384 mesh_to_volume_internal::OffsetValues<TreeType>(nodes, -offset));
3387 mesh_to_volume_internal::Renormalize<TreeType>(
3388 distTree, nodes,
buffer.get(), voxelSize));
3391 mesh_to_volume_internal::MinCombine<TreeType>(nodes,
buffer.get()));
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);
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>
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]};
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];
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>
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);
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;
3758 IntTreeT mLastPrimTree;
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) {
3904 nijk = ijk + util::COORD_OFFSETS[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(
4005 if (!(
std::abs(divisor) > 0.0))
return false;
4009 double inv_divisor = 1.0 /
divisor;
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);
4243 #endif // OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
Vec2< T > minComponent(const Vec2< T > &v1, const Vec2< T > &v2)
Return component-wise minimum of the two vectors.
void parallel_for(int64_t start, int64_t end, std::function< void(int64_t index)> &&task, parallel_options opt=parallel_options(0, Split_Y, 1))
GU_API exint quads(GU_Detail *gdp, UT_Array< GA_OffsetArray > &rings, UT_Array< GA_OffsetArray > &originalRings, GA_PrimitiveGroup *patchgroup, GA_PrimitiveGroup *loopgroup, bool smooth, fpreal smoothstrength, bool edgeloop, fpreal edgeloopPercentage)
GA_API const UT_StringHolder dist
GLboolean GLboolean GLboolean b
void clear() override
Remove all nodes from this cache, then reinsert the root node.
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Type Pow2(Type x)
Return x2.
NodeType * getNode()
Return the cached node of type NodeType. [Mainly for internal use].
IMF_EXPORT IMATH_NAMESPACE::V3f direction(const IMATH_NAMESPACE::Box2i &dataWindow, const IMATH_NAMESPACE::V2f &pixelPosition)
GLuint GLsizei const GLuint const GLintptr * offsets
OPENVDB_API const Coord COORD_OFFSETS[26]
coordinate offset table for neighboring voxels
vfloat4 sqrt(const vfloat4 &a)
#define OPENVDB_USE_VERSION_NAMESPACE
#define OPENVDB_LOG_DEBUG(mesg)
Dummy NOOP interrupter class defining interface.
T dot(const Vec3< T > &v) const
Dot product.
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...
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat s1
GLubyte GLubyte GLubyte GLubyte w
void clear()
Remove all tiles from this tree and all nodes other than the root node.
GLuint GLenum GLenum transform
SYS_FORCE_INLINE const_iterator end() const
void OIIO_API split(string_view str, std::vector< string_view > &result, string_view sep=string_view(), int maxsplit=-1)
GLint GLint GLsizei GLint GLenum GLenum type
ImageBuf OIIO_API min(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
IMATH_INTERNAL_NAMESPACE_HEADER_ENTER T abs(T a)
math::Vec4< Index32 > Vec4I
const Vec3T & min() const
Return a const reference to the minimum point of this bounding box.
const ValueT & getValue() const
Return the tile or voxel value to which this iterator is currently pointing.
GLsizei GLsizei GLfloat distance
void merge(Tree &other, MergePolicy=MERGE_ACTIVE_STATES)
Efficiently merge another tree into this tree using one of several schemes.
LeafIter beginLeaf()
Return an iterator over all leaf nodes in this tree.
fpreal64 dot(const CE_VectorT< T > &a, const CE_VectorT< T > &b)
GLboolean GLboolean GLboolean GLboolean a
Efficient multi-threaded replacement of the background values in tree.
float Sqrt(float x)
Return the square root of a floating-point value.
ImageBuf OIIO_API max(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
Defined various multi-threaded utility functions for trees.
GLdouble GLdouble GLdouble z
Tree< typename RootNodeType::template ValueConverter< Int32 >::Type > Type
bool operator<(const GU_TetrahedronFacet &a, const GU_TetrahedronFacet &b)
Vec3< T > cross(const Vec3< T > &v) const
Return the cross product of "this" vector and v;.
typedef int(WINAPI *PFNWGLRELEASEPBUFFERDCARBPROC)(HPBUFFERARB hPbuffer
OPENVDB_API const Index32 INVALID_IDX
LeafNodeT * probeLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains voxel (x, y, z), or nullptr if no such node exists...
const Vec3T & max() const
Return a const reference to the maximum point of this bounding box.
Propagate the signs of distance values from the active voxels in the narrow band to the inactive valu...
GLuint GLdouble GLdouble GLint GLint const GLdouble * points
GLuint GLsizei GLsizei * length
Axis-aligned bounding box.
typename RootNodeType::LeafNodeType LeafNodeType
GA_API const UT_StringHolder up
bool cancelGroupExecution()
math::Vec3< Index32 > Vec3I
Base class for tree-traversal iterators over tile and voxel values.
Vec2< T > maxComponent(const Vec2< T > &v1, const Vec2< T > &v2)
Return component-wise maximum of the two vectors.
GA_API const UT_StringHolder N
_RootNodeType RootNodeType
GLsizei const GLfloat * value
LeafNodeType * probeLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains voxel (x, y, z). If no such node exists, return nullptr.
Convert polygonal meshes that consist of quads and/or triangles into signed or unsigned distance fiel...
bool wasInterrupted(T *i, int percent=-1)
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
void sort(I begin, I end, const Pred &pred)
bool probeValue(const Coord &xyz, ValueType &value) const
Return the active state of the voxel as well as its value.
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
#define OPENVDB_THROW(exception, message)
void clearAllAccessors()
Clear all registered accessors.
Base class for tree-traversal iterators over all leaf nodes (but not leaf voxels) ...
Real GodunovsNormSqrd(bool isOutside, Real dP_xm, Real dP_xp, Real dP_ym, Real dP_yp, Real dP_zm, Real dP_zp)