16 #ifndef OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
17 #define OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
33 #include <tbb/blocked_range.h>
34 #include <tbb/enumerable_thread_specific.h>
35 #include <tbb/parallel_for.h>
36 #include <tbb/parallel_reduce.h>
37 #include <tbb/partitioner.h>
38 #include <tbb/task_group.h>
39 #include <tbb/task_arena.h>
47 #include <type_traits>
128 template <
typename Gr
idType,
typename MeshDataAdapter,
typename InteriorTest = std::
nullptr_t>
129 typename GridType::Ptr
133 float exteriorBandWidth = 3.0
f,
134 float interiorBandWidth = 3.0
f,
137 InteriorTest interiorTest =
nullptr,
159 template <
typename Gr
idType,
typename MeshDataAdapter,
typename Interrupter,
typename InteriorTest = std::
nullptr_t>
160 typename GridType::Ptr
162 Interrupter& interrupter,
165 float exteriorBandWidth = 3.0
f,
166 float interiorBandWidth = 3.0
f,
169 InteriorTest interiorTest =
nullptr,
186 template<
typename Po
intType,
typename PolygonType>
190 const std::vector<PolygonType>& polygons)
191 : mPointArray(points.empty() ? nullptr : &points[0])
192 , mPointArraySize(points.
size())
193 , mPolygonArray(polygons.empty() ? nullptr : &polygons[0])
194 , mPolygonArraySize(polygons.
size())
199 const PolygonType* polygonArray,
size_t polygonArraySize)
200 : mPointArray(pointArray)
201 , mPointArraySize(pointArraySize)
202 , mPolygonArray(polygonArray)
203 , mPolygonArraySize(polygonArraySize)
218 const PointType& p = mPointArray[mPolygonArray[
n][
int(v)]];
219 pos[0] = double(p[0]);
220 pos[1] = double(p[1]);
221 pos[2] = double(p[2]);
225 PointType
const *
const mPointArray;
226 size_t const mPointArraySize;
227 PolygonType
const *
const mPolygonArray;
228 size_t const mPolygonArraySize;
256 template<
typename Gr
idType>
257 typename GridType::Ptr
259 const openvdb::math::Transform& xform,
260 const std::vector<Vec3s>&
points,
261 const std::vector<Vec3I>& triangles,
262 float halfWidth =
float(LEVEL_SET_HALF_WIDTH));
265 template<
typename Gr
idType,
typename Interrupter>
266 typename GridType::Ptr
268 Interrupter& interrupter,
269 const openvdb::math::Transform& xform,
270 const std::vector<Vec3s>& points,
271 const std::vector<Vec3I>& triangles,
272 float halfWidth =
float(LEVEL_SET_HALF_WIDTH));
290 template<
typename Gr
idType>
291 typename GridType::Ptr
293 const openvdb::math::Transform& xform,
294 const std::vector<Vec3s>& points,
295 const std::vector<Vec4I>&
quads,
296 float halfWidth =
float(LEVEL_SET_HALF_WIDTH));
299 template<
typename Gr
idType,
typename Interrupter>
300 typename GridType::Ptr
302 Interrupter& interrupter,
303 const openvdb::math::Transform& xform,
304 const std::vector<Vec3s>& points,
305 const std::vector<Vec4I>& quads,
306 float halfWidth =
float(LEVEL_SET_HALF_WIDTH));
325 template<
typename Gr
idType>
326 typename GridType::Ptr
328 const openvdb::math::Transform& xform,
329 const std::vector<Vec3s>& points,
330 const std::vector<Vec3I>& triangles,
331 const std::vector<Vec4I>& quads,
332 float halfWidth =
float(LEVEL_SET_HALF_WIDTH));
335 template<
typename Gr
idType,
typename Interrupter>
336 typename GridType::Ptr
338 Interrupter& interrupter,
339 const openvdb::math::Transform& xform,
340 const std::vector<Vec3s>& points,
341 const std::vector<Vec3I>& triangles,
342 const std::vector<Vec4I>& quads,
343 float halfWidth =
float(LEVEL_SET_HALF_WIDTH));
364 template<
typename Gr
idType>
365 typename GridType::Ptr
367 const openvdb::math::Transform& xform,
368 const std::vector<Vec3s>& points,
369 const std::vector<Vec3I>& triangles,
370 const std::vector<Vec4I>& quads,
375 template<
typename Gr
idType,
typename Interrupter>
376 typename GridType::Ptr
378 Interrupter& interrupter,
379 const openvdb::math::Transform& xform,
380 const std::vector<Vec3s>& points,
381 const std::vector<Vec3I>& triangles,
382 const std::vector<Vec4I>& quads,
401 template<
typename Gr
idType>
402 typename GridType::Ptr
404 const openvdb::math::Transform& xform,
405 const std::vector<Vec3s>& points,
406 const std::vector<Vec3I>& triangles,
407 const std::vector<Vec4I>& quads,
411 template<
typename Gr
idType,
typename Interrupter>
412 typename GridType::Ptr
414 Interrupter& interrupter,
415 const openvdb::math::Transform& xform,
416 const std::vector<Vec3s>& points,
417 const std::vector<Vec3I>& triangles,
418 const std::vector<Vec4I>& quads,
431 template<
typename Gr
idType,
typename VecType>
432 typename GridType::Ptr
434 const openvdb::math::Transform& xform,
447 template <
typename FloatTreeT>
508 void convert(
const std::vector<Vec3s>& pointList,
const std::vector<Vec4I>& polygonList);
514 std::vector<Vec3d>& points, std::vector<Index32>& primitives);
535 namespace mesh_to_volume_internal {
537 template<
typename Po
intType>
538 struct TransformPoints {
540 TransformPoints(
const PointType* pointsIn, PointType* pointsOut,
541 const math::Transform& xform)
542 : mPointsIn(pointsIn), mPointsOut(pointsOut), mXform(&xform)
546 void operator()(
const tbb::blocked_range<size_t>&
range)
const {
550 for (
size_t n = range.begin(),
N = range.
end();
n <
N; ++
n) {
552 const PointType& wsP = mPointsIn[
n];
553 pos[0] = double(wsP[0]);
554 pos[1] = double(wsP[1]);
555 pos[2] = double(wsP[2]);
557 pos = mXform->worldToIndex(pos);
559 PointType& isP = mPointsOut[
n];
566 PointType
const *
const mPointsIn;
567 PointType *
const mPointsOut;
568 math::Transform
const *
const mXform;
572 template<
typename ValueType>
583 template<
typename TreeType>
584 class CombineLeafNodes
590 using LeafNodeType =
typename TreeType::LeafNodeType;
591 using Int32LeafNodeType =
typename Int32TreeType::LeafNodeType;
593 CombineLeafNodes(TreeType& lhsDistTree, Int32TreeType& lhsIdxTree,
594 LeafNodeType ** rhsDistNodes, Int32LeafNodeType ** rhsIdxNodes)
595 : mDistTree(&lhsDistTree)
596 , mIdxTree(&lhsIdxTree)
597 , mRhsDistNodes(rhsDistNodes)
598 , mRhsIdxNodes(rhsIdxNodes)
602 void operator()(
const tbb::blocked_range<size_t>& range)
const {
604 tree::ValueAccessor<TreeType> distAcc(*mDistTree);
605 tree::ValueAccessor<Int32TreeType> idxAcc(*mIdxTree);
610 for (
size_t n = range.begin(),
N = range.
end();
n <
N; ++
n) {
612 const Coord& origin = mRhsDistNodes[
n]->origin();
614 LeafNodeType* lhsDistNode = distAcc.probeLeaf(origin);
615 Int32LeafNodeType* lhsIdxNode = idxAcc.probeLeaf(origin);
617 DistValueType* lhsDistData = lhsDistNode->buffer().data();
618 IndexValueType* lhsIdxData = lhsIdxNode->buffer().data();
620 const DistValueType* rhsDistData = mRhsDistNodes[
n]->buffer().data();
621 const IndexValueType* rhsIdxData = mRhsIdxNodes[
n]->buffer().data();
628 const DistValueType& lhsValue = lhsDistData[
offset];
629 const DistValueType& rhsValue = rhsDistData[
offset];
631 if (rhsValue < lhsValue) {
632 lhsDistNode->setValueOn(
offset, rhsValue);
635 lhsIdxNode->setValueOn(
offset,
641 delete mRhsDistNodes[
n];
642 delete mRhsIdxNodes[
n];
648 TreeType *
const mDistTree;
649 Int32TreeType *
const mIdxTree;
651 LeafNodeType **
const mRhsDistNodes;
652 Int32LeafNodeType **
const mRhsIdxNodes;
659 template<
typename TreeType>
660 struct StashOriginAndStoreOffset
662 using LeafNodeType =
typename TreeType::LeafNodeType;
664 StashOriginAndStoreOffset(std::vector<LeafNodeType*>& nodes, Coord* coordinates)
665 : mNodes(nodes.empty() ? nullptr : &nodes[0]), mCoordinates(coordinates)
669 void operator()(
const tbb::blocked_range<size_t>& range)
const {
670 for (
size_t n = range.begin(), N = range.end();
n <
N; ++
n) {
671 Coord& origin =
const_cast<Coord&
>(mNodes[
n]->origin());
672 mCoordinates[
n] = origin;
673 origin[0] =
static_cast<int>(
n);
677 LeafNodeType **
const mNodes;
678 Coord *
const mCoordinates;
682 template<
typename TreeType>
685 using LeafNodeType =
typename TreeType::LeafNodeType;
687 RestoreOrigin(std::vector<LeafNodeType*>& nodes,
const Coord* coordinates)
688 : mNodes(nodes.empty() ? nullptr : &nodes[0]), mCoordinates(coordinates)
692 void operator()(
const tbb::blocked_range<size_t>& range)
const {
693 for (
size_t n = range.begin(), N = range.end();
n <
N; ++
n) {
694 Coord& origin =
const_cast<Coord&
>(mNodes[
n]->origin());
695 origin[0] = mCoordinates[
n][0];
699 LeafNodeType **
const mNodes;
700 Coord
const *
const mCoordinates;
704 template<
typename TreeType>
705 class ComputeNodeConnectivity
708 using LeafNodeType =
typename TreeType::LeafNodeType;
710 ComputeNodeConnectivity(
const TreeType& tree,
const Coord* coordinates,
713 , mCoordinates(coordinates)
715 , mNumNodes(numNodes)
720 ComputeNodeConnectivity(
const ComputeNodeConnectivity&) =
default;
723 ComputeNodeConnectivity&
operator=(
const ComputeNodeConnectivity&) =
delete;
725 void operator()(
const tbb::blocked_range<size_t>& range)
const {
727 size_t* offsetsNextX = mOffsets;
728 size_t* offsetsPrevX = mOffsets + mNumNodes;
729 size_t* offsetsNextY = mOffsets + mNumNodes * 2;
730 size_t* offsetsPrevY = mOffsets + mNumNodes * 3;
731 size_t* offsetsNextZ = mOffsets + mNumNodes * 4;
732 size_t* offsetsPrevZ = mOffsets + mNumNodes * 5;
734 tree::ValueAccessor<const TreeType> acc(*mTree);
736 const Int32 DIM =
static_cast<Int32>(LeafNodeType::DIM);
738 for (
size_t n = range.begin(), N = range.end();
n <
N; ++
n) {
739 const Coord& origin = mCoordinates[
n];
740 offsetsNextX[
n] = findNeighbourNode(acc, origin, Coord(DIM, 0, 0));
741 offsetsPrevX[
n] = findNeighbourNode(acc, origin, Coord(-DIM, 0, 0));
742 offsetsNextY[
n] = findNeighbourNode(acc, origin, Coord(0, DIM, 0));
743 offsetsPrevY[
n] = findNeighbourNode(acc, origin, Coord(0, -DIM, 0));
744 offsetsNextZ[
n] = findNeighbourNode(acc, origin, Coord(0, 0, DIM));
745 offsetsPrevZ[
n] = findNeighbourNode(acc, origin, Coord(0, 0, -DIM));
749 size_t findNeighbourNode(tree::ValueAccessor<const TreeType>& acc,
750 const Coord&
start,
const Coord& step)
const
752 Coord ijk = start + step;
755 while (bbox.isInside(ijk)) {
756 const LeafNodeType* node = acc.probeConstLeaf(ijk);
757 if (node)
return static_cast<size_t>(node->origin()[0]);
766 TreeType
const *
const mTree;
767 Coord
const *
const mCoordinates;
768 size_t *
const mOffsets;
770 const size_t mNumNodes;
775 template<
typename TreeType>
776 struct LeafNodeConnectivityTable
780 using LeafNodeType =
typename TreeType::LeafNodeType;
782 LeafNodeConnectivityTable(TreeType& tree)
784 mLeafNodes.reserve(tree.leafCount());
785 tree.getNodes(mLeafNodes);
787 if (mLeafNodes.empty())
return;
790 tree.evalLeafBoundingBox(bbox);
792 const tbb::blocked_range<size_t>
range(0, mLeafNodes.size());
796 std::unique_ptr<Coord[]> coordinates{
new Coord[mLeafNodes.size()]};
798 StashOriginAndStoreOffset<TreeType>(mLeafNodes, coordinates.get()));
801 mOffsets.reset(
new size_t[mLeafNodes.size() * 6]);
805 tree, coordinates.get(), mOffsets.get(), mLeafNodes.size(), bbox));
808 tbb::parallel_for(range, RestoreOrigin<TreeType>(mLeafNodes, coordinates.get()));
811 size_t size()
const {
return mLeafNodes.size(); }
813 std::vector<LeafNodeType*>& nodes() {
return mLeafNodes; }
814 const std::vector<LeafNodeType*>& nodes()
const {
return mLeafNodes; }
817 const size_t* offsetsNextX()
const {
return mOffsets.get(); }
818 const size_t* offsetsPrevX()
const {
return mOffsets.get() + mLeafNodes.size(); }
820 const size_t* offsetsNextY()
const {
return mOffsets.get() + mLeafNodes.size() * 2; }
821 const size_t* offsetsPrevY()
const {
return mOffsets.get() + mLeafNodes.size() * 3; }
823 const size_t* offsetsNextZ()
const {
return mOffsets.get() + mLeafNodes.size() * 4; }
824 const size_t* offsetsPrevZ()
const {
return mOffsets.get() + mLeafNodes.size() * 5; }
827 std::vector<LeafNodeType*> mLeafNodes;
828 std::unique_ptr<size_t[]> mOffsets;
832 template<
typename TreeType>
833 class SweepExteriorSign
840 using LeafNodeType =
typename TreeType::LeafNodeType;
841 using ConnectivityTable = LeafNodeConnectivityTable<TreeType>;
843 SweepExteriorSign(
Axis axis,
const std::vector<size_t>& startNodeIndices,
844 ConnectivityTable& connectivity)
845 : mStartNodeIndices(startNodeIndices.empty() ? nullptr : &startNodeIndices[0])
846 , mConnectivity(&connectivity)
851 void operator()(
const tbb::blocked_range<size_t>& range)
const {
853 constexpr
Int32 DIM =
static_cast<Int32>(LeafNodeType::DIM);
855 std::vector<LeafNodeType*>& nodes = mConnectivity->nodes();
858 size_t idxA = 0, idxB = 1;
861 const size_t* nextOffsets = mConnectivity->offsetsNextZ();
862 const size_t* prevOffsets = mConnectivity->offsetsPrevZ();
870 nextOffsets = mConnectivity->offsetsNextY();
871 prevOffsets = mConnectivity->offsetsPrevY();
873 }
else if (mAxis ==
X_AXIS) {
879 nextOffsets = mConnectivity->offsetsNextX();
880 prevOffsets = mConnectivity->offsetsPrevX();
888 for (
size_t n = range.begin(), N = range.end();
n <
N; ++
n) {
890 size_t startOffset = mStartNodeIndices[
n];
895 for (a = 0; a < DIM; ++
a) {
896 for (b = 0; b < DIM; ++
b) {
898 pos =
static_cast<Int32>(LeafNodeType::coordToOffset(ijk));
899 size_t offset = startOffset;
902 while ( offset != ConnectivityTable::INVALID_OFFSET &&
903 traceVoxelLine(*nodes[offset], pos, step) ) {
906 offset = nextOffsets[
offset];
911 while (offset != ConnectivityTable::INVALID_OFFSET) {
913 offset = nextOffsets[
offset];
918 pos += step * (DIM - 1);
919 while ( offset != ConnectivityTable::INVALID_OFFSET &&
920 traceVoxelLine(*nodes[offset], pos, -step)) {
921 offset = prevOffsets[
offset];
929 bool traceVoxelLine(LeafNodeType& node, Int32 pos,
const Int32 step)
const {
933 bool isOutside =
true;
935 for (Index i = 0; i < LeafNodeType::DIM; ++i) {
944 if (!(dist >
ValueType(0.75))) isOutside =
false;
957 size_t const *
const mStartNodeIndices;
958 ConnectivityTable *
const mConnectivity;
964 template<
typename LeafNodeType>
966 seedFill(LeafNodeType& node)
969 using Queue = std::deque<Index>;
977 if (data[pos] < 0.0) seedPoints.push_back(pos);
980 if (seedPoints.empty())
return;
983 for (Queue::iterator it = seedPoints.begin(); it != seedPoints.end(); ++it) {
991 Index pos(0), nextPos(0);
993 while (!seedPoints.empty()) {
995 pos = seedPoints.back();
996 seedPoints.pop_back();
1004 ijk = LeafNodeType::offsetToLocalCoord(pos);
1007 nextPos = pos - LeafNodeType::DIM * LeafNodeType::DIM;
1008 if (data[nextPos] >
ValueType(0.75)) seedPoints.push_back(nextPos);
1011 if (ijk[0] != (LeafNodeType::DIM - 1)) {
1012 nextPos = pos + LeafNodeType::DIM * LeafNodeType::DIM;
1013 if (data[nextPos] >
ValueType(0.75)) seedPoints.push_back(nextPos);
1017 nextPos = pos - LeafNodeType::DIM;
1018 if (data[nextPos] >
ValueType(0.75)) seedPoints.push_back(nextPos);
1021 if (ijk[1] != (LeafNodeType::DIM - 1)) {
1022 nextPos = pos + LeafNodeType::DIM;
1023 if (data[nextPos] >
ValueType(0.75)) seedPoints.push_back(nextPos);
1028 if (data[nextPos] >
ValueType(0.75)) seedPoints.push_back(nextPos);
1031 if (ijk[2] != (LeafNodeType::DIM - 1)) {
1033 if (data[nextPos] >
ValueType(0.75)) seedPoints.push_back(nextPos);
1040 template<
typename LeafNodeType>
1042 scanFill(LeafNodeType& node)
1044 bool updatedNode =
false;
1051 bool updatedSign =
true;
1052 while (updatedSign) {
1054 updatedSign =
false;
1062 ijk = LeafNodeType::offsetToLocalCoord(pos);
1065 if (ijk[2] != 0 && data[pos - 1] <
ValueType(0.0)) {
1070 }
else if (ijk[2] != (LeafNodeType::DIM - 1) && data[pos + 1] <
ValueType(0.0)) {
1075 }
else if (ijk[1] != 0 && data[pos - LeafNodeType::DIM] <
ValueType(0.0)) {
1080 }
else if (ijk[1] != (LeafNodeType::DIM - 1)
1081 && data[pos + LeafNodeType::DIM] <
ValueType(0.0))
1087 }
else if (ijk[0] != 0
1088 && data[pos - LeafNodeType::DIM * LeafNodeType::DIM] <
ValueType(0.0))
1094 }
else if (ijk[0] != (LeafNodeType::DIM - 1)
1095 && data[pos + LeafNodeType::DIM * LeafNodeType::DIM] <
ValueType(0.0))
1103 updatedNode |= updatedSign;
1110 template<
typename TreeType>
1111 class SeedFillExteriorSign
1115 using LeafNodeType =
typename TreeType::LeafNodeType;
1117 SeedFillExteriorSign(std::vector<LeafNodeType*>& nodes,
const bool* changedNodeMask)
1118 : mNodes(nodes.empty() ? nullptr : &nodes[0])
1119 , mChangedNodeMask(changedNodeMask)
1123 void operator()(
const tbb::blocked_range<size_t>& range)
const {
1124 for (
size_t n = range.begin(), N = range.end();
n <
N; ++
n) {
1125 if (mChangedNodeMask[
n]) {
1131 scanFill(*mNodes[n]);
1136 LeafNodeType **
const mNodes;
1137 const bool *
const mChangedNodeMask;
1141 template<
typename ValueType>
1146 void operator()(
const tbb::blocked_range<size_t>& range)
const {
1148 for (
size_t n = range.begin(), N = range.end();
n <
N; ++
n) {
1158 template<
typename ValueType>
1162 const auto grainSize = std::max<size_t>(
1163 length / tbb::this_task_arena::max_concurrency(), 1024);
1164 const tbb::blocked_range<size_t>
range(0, length, grainSize);
1165 tbb::parallel_for(range, FillArray<ValueType>(array, val), tbb::simple_partitioner());
1169 template<
typename TreeType>
1174 using LeafNodeType =
typename TreeType::LeafNodeType;
1176 SyncVoxelMask(std::vector<LeafNodeType*>& nodes,
1177 const bool* changedNodeMask,
bool* changedVoxelMask)
1178 : mNodes(nodes.empty() ? nullptr : &nodes[0])
1179 , mChangedNodeMask(changedNodeMask)
1180 , mChangedVoxelMask(changedVoxelMask)
1184 void operator()(
const tbb::blocked_range<size_t>& range)
const {
1185 for (
size_t n = range.begin(), N = range.end();
n <
N; ++
n) {
1187 if (mChangedNodeMask[
n]) {
1190 ValueType* data = mNodes[
n]->buffer().data();
1202 LeafNodeType **
const mNodes;
1203 bool const *
const mChangedNodeMask;
1204 bool *
const mChangedVoxelMask;
1208 template<
typename TreeType>
1213 using LeafNodeType =
typename TreeType::LeafNodeType;
1214 using ConnectivityTable = LeafNodeConnectivityTable<TreeType>;
1216 SeedPoints(ConnectivityTable& connectivity,
1217 bool* changedNodeMask,
bool* nodeMask,
bool* changedVoxelMask)
1218 : mConnectivity(&connectivity)
1219 , mChangedNodeMask(changedNodeMask)
1220 , mNodeMask(nodeMask)
1221 , mChangedVoxelMask(changedVoxelMask)
1225 void operator()(
const tbb::blocked_range<size_t>& range)
const {
1227 for (
size_t n = range.begin(), N = range.end(); n <
N; ++
n) {
1228 bool changedValue =
false;
1230 changedValue |= processZ(n,
true);
1231 changedValue |= processZ(n,
false);
1233 changedValue |= processY(n,
true);
1234 changedValue |= processY(n,
false);
1236 changedValue |= processX(n,
true);
1237 changedValue |= processX(n,
false);
1239 mNodeMask[
n] = changedValue;
1244 bool processZ(
const size_t n,
bool firstFace)
const
1246 const size_t offset =
1247 firstFace ? mConnectivity->offsetsPrevZ()[
n] : mConnectivity->offsetsNextZ()[
n];
1248 if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1252 const ValueType* lhsData = mConnectivity->nodes()[
n]->buffer().data();
1253 const ValueType* rhsData = mConnectivity->nodes()[
offset]->buffer().data();
1255 const Index lastOffset = LeafNodeType::DIM - 1;
1256 const Index lhsOffset =
1257 firstFace ? 0 :
lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1259 Index tmpPos(0), pos(0);
1260 bool changedValue =
false;
1262 for (Index
x = 0;
x < LeafNodeType::DIM; ++
x) {
1263 tmpPos =
x << (2 * LeafNodeType::LOG2DIM);
1264 for (Index
y = 0;
y < LeafNodeType::DIM; ++
y) {
1265 pos = tmpPos + (
y << LeafNodeType::LOG2DIM);
1267 if (lhsData[pos + lhsOffset] >
ValueType(0.75)) {
1268 if (rhsData[pos + rhsOffset] <
ValueType(0.0)) {
1269 changedValue =
true;
1270 mask[pos + lhsOffset] =
true;
1276 return changedValue;
1282 bool processY(
const size_t n,
bool firstFace)
const
1284 const size_t offset =
1285 firstFace ? mConnectivity->offsetsPrevY()[
n] : mConnectivity->offsetsNextY()[
n];
1286 if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1290 const ValueType* lhsData = mConnectivity->nodes()[
n]->buffer().data();
1291 const ValueType* rhsData = mConnectivity->nodes()[
offset]->buffer().data();
1293 const Index lastOffset = LeafNodeType::DIM * (LeafNodeType::DIM - 1);
1294 const Index lhsOffset =
1295 firstFace ? 0 :
lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1297 Index tmpPos(0), pos(0);
1298 bool changedValue =
false;
1300 for (Index
x = 0;
x < LeafNodeType::DIM; ++
x) {
1301 tmpPos =
x << (2 * LeafNodeType::LOG2DIM);
1302 for (Index
z = 0;
z < LeafNodeType::DIM; ++
z) {
1305 if (lhsData[pos + lhsOffset] >
ValueType(0.75)) {
1306 if (rhsData[pos + rhsOffset] <
ValueType(0.0)) {
1307 changedValue =
true;
1308 mask[pos + lhsOffset] =
true;
1314 return changedValue;
1320 bool processX(
const size_t n,
bool firstFace)
const
1322 const size_t offset =
1323 firstFace ? mConnectivity->offsetsPrevX()[
n] : mConnectivity->offsetsNextX()[
n];
1324 if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1328 const ValueType* lhsData = mConnectivity->nodes()[
n]->buffer().data();
1329 const ValueType* rhsData = mConnectivity->nodes()[
offset]->buffer().data();
1331 const Index lastOffset = LeafNodeType::DIM * LeafNodeType::DIM * (LeafNodeType::DIM-1);
1332 const Index lhsOffset =
1333 firstFace ? 0 :
lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1335 Index tmpPos(0), pos(0);
1336 bool changedValue =
false;
1338 for (Index
y = 0;
y < LeafNodeType::DIM; ++
y) {
1339 tmpPos =
y << LeafNodeType::LOG2DIM;
1340 for (Index z = 0; z < LeafNodeType::DIM; ++
z) {
1343 if (lhsData[pos + lhsOffset] >
ValueType(0.75)) {
1344 if (rhsData[pos + rhsOffset] <
ValueType(0.0)) {
1345 changedValue =
true;
1346 mask[pos + lhsOffset] =
true;
1352 return changedValue;
1358 ConnectivityTable *
const mConnectivity;
1359 bool *
const mChangedNodeMask;
1360 bool *
const mNodeMask;
1361 bool *
const mChangedVoxelMask;
1367 template<
typename TreeType,
typename MeshDataAdapter>
1368 struct ComputeIntersectingVoxelSign
1371 using LeafNodeType =
typename TreeType::LeafNodeType;
1373 using Int32LeafNodeType =
typename Int32TreeType::LeafNodeType;
1376 using MaskArray = std::unique_ptr<bool[]>;
1377 using LocalData = std::pair<PointArray, MaskArray>;
1378 using LocalDataTable = tbb::enumerable_thread_specific<LocalData>;
1380 ComputeIntersectingVoxelSign(
1381 std::vector<LeafNodeType*>& distNodes,
1382 const TreeType& distTree,
1383 const Int32TreeType& indexTree,
1385 : mDistNodes(distNodes.empty() ? nullptr : &distNodes[0])
1386 , mDistTree(&distTree)
1387 , mIndexTree(&indexTree)
1389 , mLocalDataTable(new LocalDataTable())
1394 void operator()(
const tbb::blocked_range<size_t>& range)
const {
1396 tree::ValueAccessor<const TreeType> distAcc(*mDistTree);
1397 tree::ValueAccessor<const Int32TreeType> idxAcc(*mIndexTree);
1401 Index xPos(0), yPos(0);
1402 Coord ijk, nijk, nodeMin, nodeMax;
1403 Vec3d cp, xyz, nxyz, dir1, dir2;
1405 LocalData& localData = mLocalDataTable->local();
1408 if (!points) points.reset(
new Vec3d[LeafNodeType::SIZE * 2]);
1410 MaskArray& mask = localData.second;
1411 if (!mask) mask.reset(
new bool[LeafNodeType::SIZE]);
1414 typename LeafNodeType::ValueOnCIter it;
1416 for (
size_t n = range.begin(), N = range.end(); n <
N; ++
n) {
1418 LeafNodeType& node = *mDistNodes[
n];
1421 const Int32LeafNodeType* idxNode = idxAcc.probeConstLeaf(node.origin());
1422 const Int32* idxData = idxNode->buffer().data();
1424 nodeMin = node.origin();
1425 nodeMax = nodeMin.offsetBy(LeafNodeType::DIM - 1);
1428 memset(mask.get(), 0,
sizeof(
bool) * LeafNodeType::SIZE);
1430 for (it = node.cbeginValueOn(); it; ++it) {
1431 Index pos = it.pos();
1434 if (dist < 0.0 || dist > 0.75)
continue;
1436 ijk = node.offsetToGlobalCoord(pos);
1438 xyz[0] = double(ijk[0]);
1439 xyz[1] = double(ijk[1]);
1440 xyz[2] = double(ijk[2]);
1446 bool flipSign =
false;
1448 for (nijk[0] = bbox.min()[0]; nijk[0] <= bbox.max()[0] && !flipSign; ++nijk[0]) {
1449 xPos = (nijk[0] & (LeafNodeType::DIM - 1u)) << (2 * LeafNodeType::LOG2DIM);
1450 for (nijk[1]=bbox.min()[1]; nijk[1] <= bbox.max()[1] && !flipSign; ++nijk[1]) {
1451 yPos = xPos + ((nijk[1] & (LeafNodeType::DIM-1u)) << LeafNodeType::LOG2DIM);
1452 for (nijk[2] = bbox.min()[2]; nijk[2] <= bbox.max()[2]; ++nijk[2]) {
1453 pos = yPos + (nijk[2] & (LeafNodeType::DIM - 1u));
1455 const Int32& polyIdx = idxData[pos];
1460 const Index pointIndex = pos * 2;
1466 nxyz[0] = double(nijk[0]);
1467 nxyz[1] = double(nijk[1]);
1468 nxyz[2] = double(nijk[2]);
1470 Vec3d& point = points[pointIndex];
1472 point = closestPoint(nxyz, polyIdx);
1475 direction = nxyz - point;
1476 direction.normalize();
1479 dir1 = xyz - points[pointIndex];
1482 if (points[pointIndex + 1].
dot(dir1) > 0.0) {
1493 for (Int32 m = 0; m < 26; ++m) {
1496 if (!bbox.isInside(nijk) && distAcc.probeValue(nijk, nval) && nval<-0.75) {
1497 nxyz[0] = double(nijk[0]);
1498 nxyz[1] = double(nijk[1]);
1499 nxyz[2] = double(nijk[2]);
1501 cp = closestPoint(nxyz, idxAcc.getValue(nijk));
1509 if (dir2.dot(dir1) > 0.0) {
1523 Vec3d closestPoint(
const Vec3d& center, Int32 polyIdx)
const
1527 const size_t polygon = size_t(polyIdx);
1528 mMesh->getIndexSpacePoint(polygon, 0, a);
1529 mMesh->getIndexSpacePoint(polygon, 1, b);
1530 mMesh->getIndexSpacePoint(polygon, 2, c);
1534 if (4 == mMesh->vertexCount(polygon)) {
1536 mMesh->getIndexSpacePoint(polygon, 3, b);
1540 if ((center - c).lengthSqr() < (center - cp).lengthSqr()) {
1549 LeafNodeType **
const mDistNodes;
1550 TreeType
const *
const mDistTree;
1551 Int32TreeType
const *
const mIndexTree;
1554 SharedPtr<LocalDataTable> mLocalDataTable;
1561 template<
typename LeafNodeType>
1563 maskNodeInternalNeighbours(
const Index pos,
bool (&mask)[26])
1565 using NodeT = LeafNodeType;
1567 const Coord ijk = NodeT::offsetToLocalCoord(pos);
1571 mask[0] = ijk[0] != (NodeT::DIM - 1);
1573 mask[1] = ijk[0] != 0;
1575 mask[2] = ijk[1] != (NodeT::DIM - 1);
1577 mask[3] = ijk[1] != 0;
1579 mask[4] = ijk[2] != (NodeT::DIM - 1);
1581 mask[5] = ijk[2] != 0;
1585 mask[6] = mask[0] && mask[5];
1587 mask[7] = mask[1] && mask[5];
1589 mask[8] = mask[0] && mask[4];
1591 mask[9] = mask[1] && mask[4];
1593 mask[10] = mask[0] && mask[2];
1595 mask[11] = mask[1] && mask[2];
1597 mask[12] = mask[0] && mask[3];
1599 mask[13] = mask[1] && mask[3];
1601 mask[14] = mask[3] && mask[4];
1603 mask[15] = mask[3] && mask[5];
1605 mask[16] = mask[2] && mask[4];
1607 mask[17] = mask[2] && mask[5];
1611 mask[18] = mask[1] && mask[3] && mask[5];
1613 mask[19] = mask[1] && mask[3] && mask[4];
1615 mask[20] = mask[0] && mask[3] && mask[4];
1617 mask[21] = mask[0] && mask[3] && mask[5];
1619 mask[22] = mask[1] && mask[2] && mask[5];
1621 mask[23] = mask[1] && mask[2] && mask[4];
1623 mask[24] = mask[0] && mask[2] && mask[4];
1625 mask[25] = mask[0] && mask[2] && mask[5];
1629 template<
typename Compare,
typename LeafNodeType>
1633 using NodeT = LeafNodeType;
1636 if (mask[5] && Compare::check(data[pos - 1]))
return true;
1638 if (mask[4] && Compare::check(data[pos + 1]))
return true;
1640 if (mask[3] && Compare::check(data[pos - NodeT::DIM]))
return true;
1642 if (mask[2] && Compare::check(data[pos + NodeT::DIM]))
return true;
1644 if (mask[1] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM]))
return true;
1646 if (mask[0] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM]))
return true;
1648 if (mask[6] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM]))
return true;
1650 if (mask[7] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - 1]))
return true;
1652 if (mask[8] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + 1]))
return true;
1654 if (mask[9] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + 1]))
return true;
1656 if (mask[10] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM]))
return true;
1658 if (mask[11] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM]))
return true;
1660 if (mask[12] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM]))
return true;
1662 if (mask[13] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM]))
return true;
1664 if (mask[14] && Compare::check(data[pos - NodeT::DIM + 1]))
return true;
1666 if (mask[15] && Compare::check(data[pos - NodeT::DIM - 1]))
return true;
1668 if (mask[16] && Compare::check(data[pos + NodeT::DIM + 1]))
return true;
1670 if (mask[17] && Compare::check(data[pos + NodeT::DIM - 1]))
return true;
1672 if (mask[18] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM - 1]))
return true;
1674 if (mask[19] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM + 1]))
return true;
1676 if (mask[20] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM + 1]))
return true;
1678 if (mask[21] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM - 1]))
return true;
1680 if (mask[22] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM - 1]))
return true;
1682 if (mask[23] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM + 1]))
return true;
1684 if (mask[24] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM + 1]))
return true;
1686 if (mask[25] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM - 1]))
return true;
1692 template<
typename Compare,
typename AccessorType>
1694 checkNeighbours(
const Coord& ijk, AccessorType& acc,
bool (&mask)[26])
1696 for (Int32 m = 0; m < 26; ++m) {
1697 if (!mask[m] && Compare::check(acc.getValue(ijk + util::COORD_OFFSETS[m]))) {
1706 template<
typename TreeType>
1707 struct ValidateIntersectingVoxels
1710 using LeafNodeType =
typename TreeType::LeafNodeType;
1712 struct IsNegative {
static bool check(
const ValueType v) {
return v <
ValueType(0.0); } };
1714 ValidateIntersectingVoxels(TreeType& tree, std::vector<LeafNodeType*>& nodes)
1716 , mNodes(nodes.empty() ? nullptr : &nodes[0])
1720 void operator()(
const tbb::blocked_range<size_t>& range)
const
1722 tree::ValueAccessor<const TreeType> acc(*mTree);
1723 bool neighbourMask[26];
1725 for (
size_t n = range.begin(), N = range.end(); n <
N; ++
n) {
1727 LeafNodeType& node = *mNodes[
n];
1730 typename LeafNodeType::ValueOnCIter it;
1731 for (it = node.cbeginValueOn(); it; ++it) {
1733 const Index pos = it.pos();
1736 if (dist < 0.0 || dist > 0.75)
continue;
1739 maskNodeInternalNeighbours<LeafNodeType>(pos, neighbourMask);
1741 const bool hasNegativeNeighbour =
1742 checkNeighbours<IsNegative, LeafNodeType>(pos,
data, neighbourMask) ||
1743 checkNeighbours<IsNegative>(node.offsetToGlobalCoord(pos), acc, neighbourMask);
1745 if (!hasNegativeNeighbour) {
1747 dist =
ValueType(0.75) + Tolerance<ValueType>::epsilon();
1753 TreeType *
const mTree;
1754 LeafNodeType **
const mNodes;
1758 template<
typename TreeType>
1759 struct RemoveSelfIntersectingSurface
1762 using LeafNodeType =
typename TreeType::LeafNodeType;
1767 RemoveSelfIntersectingSurface(std::vector<LeafNodeType*>& nodes,
1768 TreeType& distTree, Int32TreeType& indexTree)
1769 : mNodes(nodes.empty() ? nullptr : &nodes[0])
1770 , mDistTree(&distTree)
1771 , mIndexTree(&indexTree)
1775 void operator()(
const tbb::blocked_range<size_t>& range)
const
1777 tree::ValueAccessor<const TreeType> distAcc(*mDistTree);
1778 tree::ValueAccessor<Int32TreeType> idxAcc(*mIndexTree);
1779 bool neighbourMask[26];
1781 for (
size_t n = range.begin(), N = range.end(); n <
N; ++
n) {
1783 LeafNodeType& distNode = *mNodes[
n];
1784 ValueType* data = distNode.buffer().data();
1786 typename Int32TreeType::LeafNodeType* idxNode =
1787 idxAcc.probeLeaf(distNode.origin());
1789 typename LeafNodeType::ValueOnCIter it;
1790 for (it = distNode.cbeginValueOn(); it; ++it) {
1792 const Index pos = it.pos();
1794 if (!(data[pos] > 0.75))
continue;
1797 maskNodeInternalNeighbours<LeafNodeType>(pos, neighbourMask);
1799 const bool hasBoundaryNeighbour =
1800 checkNeighbours<Comp, LeafNodeType>(pos,
data, neighbourMask) ||
1801 checkNeighbours<Comp>(distNode.offsetToGlobalCoord(pos),distAcc,neighbourMask);
1803 if (!hasBoundaryNeighbour) {
1804 distNode.setValueOff(pos);
1805 idxNode->setValueOff(pos);
1811 LeafNodeType * *
const mNodes;
1812 TreeType *
const mDistTree;
1813 Int32TreeType *
const mIndexTree;
1820 template<
typename NodeType>
1821 struct ReleaseChildNodes
1823 ReleaseChildNodes(NodeType ** nodes) : mNodes(nodes) {}
1825 void operator()(
const tbb::blocked_range<size_t>& range)
const {
1827 using NodeMaskType =
typename NodeType::NodeMaskType;
1829 for (
size_t n = range.begin(), N = range.end(); n <
N; ++
n) {
1830 const_cast<NodeMaskType&
>(mNodes[
n]->getChildMask()).setOff();
1834 NodeType **
const mNodes;
1838 template<
typename TreeType>
1840 releaseLeafNodes(TreeType& tree)
1842 using RootNodeType =
typename TreeType::RootNodeType;
1843 using NodeChainType =
typename RootNodeType::NodeChainType;
1844 using InternalNodeType =
typename NodeChainType::template Get<1>;
1846 std::vector<InternalNodeType*> nodes;
1847 tree.getNodes(nodes);
1850 ReleaseChildNodes<InternalNodeType>(nodes.empty() ?
nullptr : &nodes[0]));
1854 template<
typename TreeType>
1855 struct StealUniqueLeafNodes
1857 using LeafNodeType =
typename TreeType::LeafNodeType;
1859 StealUniqueLeafNodes(TreeType& lhsTree, TreeType& rhsTree,
1860 std::vector<LeafNodeType*>& overlappingNodes)
1861 : mLhsTree(&lhsTree)
1862 , mRhsTree(&rhsTree)
1863 , mNodes(&overlappingNodes)
1867 void operator()()
const {
1869 std::vector<LeafNodeType*> rhsLeafNodes;
1871 rhsLeafNodes.reserve(mRhsTree->leafCount());
1874 mRhsTree->stealNodes(rhsLeafNodes);
1876 tree::ValueAccessor<TreeType> acc(*mLhsTree);
1878 for (
size_t n = 0, N = rhsLeafNodes.size(); n <
N; ++
n) {
1879 if (!acc.probeLeaf(rhsLeafNodes[n]->origin())) {
1880 acc.addLeaf(rhsLeafNodes[n]);
1882 mNodes->push_back(rhsLeafNodes[n]);
1888 TreeType *
const mLhsTree;
1889 TreeType *
const mRhsTree;
1890 std::vector<LeafNodeType*> *
const mNodes;
1894 template<
typename DistTreeType,
typename IndexTreeType>
1896 combineData(DistTreeType& lhsDist, IndexTreeType& lhsIdx,
1897 DistTreeType& rhsDist, IndexTreeType& rhsIdx)
1899 using DistLeafNodeType =
typename DistTreeType::LeafNodeType;
1900 using IndexLeafNodeType =
typename IndexTreeType::LeafNodeType;
1902 std::vector<DistLeafNodeType*> overlappingDistNodes;
1903 std::vector<IndexLeafNodeType*> overlappingIdxNodes;
1906 tbb::task_group tasks;
1907 tasks.run(StealUniqueLeafNodes<DistTreeType>(lhsDist, rhsDist, overlappingDistNodes));
1908 tasks.run(StealUniqueLeafNodes<IndexTreeType>(lhsIdx, rhsIdx, overlappingIdxNodes));
1912 if (!overlappingDistNodes.empty() && !overlappingIdxNodes.empty()) {
1914 CombineLeafNodes<DistTreeType>(lhsDist, lhsIdx,
1915 &overlappingDistNodes[0], &overlappingIdxNodes[0]));
1925 template<
typename TreeType>
1926 struct VoxelizationData {
1928 using Ptr = std::unique_ptr<VoxelizationData>;
1934 using FloatTreeAcc = tree::ValueAccessor<TreeType>;
1935 using Int32TreeAcc = tree::ValueAccessor<Int32TreeType>;
1936 using UCharTreeAcc = tree::ValueAccessor<UCharTreeType>;
1943 , indexAcc(indexTree)
1944 , primIdTree(MaxPrimId)
1945 , primIdAcc(primIdTree)
1951 FloatTreeAcc distAcc;
1953 Int32TreeType indexTree;
1954 Int32TreeAcc indexAcc;
1956 UCharTreeType primIdTree;
1957 UCharTreeAcc primIdAcc;
1959 unsigned char getNewPrimId() {
1975 if (mPrimCount == MaxPrimId || primIdTree.leafCount() > 1000) {
1977 primIdTree.root().clear();
1978 primIdTree.clearAllAccessors();
1982 return mPrimCount++;
1987 enum { MaxPrimId = 100 };
1989 unsigned char mPrimCount;
1993 template<
typename TreeType,
typename MeshDataAdapter,
typename Interrupter = util::NullInterrupter>
1994 class VoxelizePolygons
1998 using VoxelizationDataType = VoxelizationData<TreeType>;
1999 using DataTable = tbb::enumerable_thread_specific<typename VoxelizationDataType::Ptr>;
2001 VoxelizePolygons(DataTable& dataTable,
2003 Interrupter* interrupter =
nullptr)
2004 : mDataTable(&dataTable)
2006 , mInterrupter(interrupter)
2010 void operator()(
const tbb::blocked_range<size_t>& range)
const {
2012 typename VoxelizationDataType::Ptr& dataPtr = mDataTable->local();
2013 if (!dataPtr) dataPtr.reset(
new VoxelizationDataType());
2017 for (
size_t n = range.begin(), N = range.end(); n <
N; ++
n) {
2020 thread::cancelGroupExecution();
2024 const size_t numVerts = mMesh->vertexCount(n);
2027 if (numVerts == 3 || numVerts == 4) {
2029 prim.index =
Int32(n);
2031 mMesh->getIndexSpacePoint(n, 0, prim.a);
2032 mMesh->getIndexSpacePoint(n, 1, prim.b);
2033 mMesh->getIndexSpacePoint(n, 2, prim.c);
2035 evalTriangle(prim, *dataPtr);
2037 if (numVerts == 4) {
2038 mMesh->getIndexSpacePoint(n, 3, prim.b);
2039 evalTriangle(prim, *dataPtr);
2047 bool wasInterrupted()
const {
return mInterrupter && mInterrupter->wasInterrupted(); }
2053 enum { POLYGON_LIMIT = 1000 };
2055 SubTask(
const Triangle& prim, DataTable& dataTable,
2056 int subdivisionCount,
size_t polygonCount,
2057 Interrupter* interrupter =
nullptr)
2058 : mLocalDataTable(&dataTable)
2060 , mSubdivisionCount(subdivisionCount)
2061 , mPolygonCount(polygonCount)
2062 , mInterrupter(interrupter)
2066 void operator()()
const
2068 if (mSubdivisionCount <= 0 || mPolygonCount >= POLYGON_LIMIT) {
2070 typename VoxelizationDataType::Ptr& dataPtr = mLocalDataTable->local();
2071 if (!dataPtr) dataPtr.reset(
new VoxelizationDataType());
2073 voxelizeTriangle(mPrim, *dataPtr, mInterrupter);
2075 }
else if (!(mInterrupter && mInterrupter->wasInterrupted())) {
2076 spawnTasks(mPrim, *mLocalDataTable, mSubdivisionCount, mPolygonCount, mInterrupter);
2080 DataTable *
const mLocalDataTable;
2081 Triangle
const mPrim;
2082 int const mSubdivisionCount;
2083 size_t const mPolygonCount;
2084 Interrupter *
const mInterrupter;
2087 inline static int evalSubdivisionCount(
const Triangle& prim)
2089 const double ax = prim.a[0], bx = prim.b[0], cx = prim.c[0];
2092 const double ay = prim.a[1], by = prim.b[1], cy = prim.c[1];
2095 const double az = prim.a[2], bz = prim.b[2], cz = prim.c[2];
2101 void evalTriangle(
const Triangle& prim, VoxelizationDataType& data)
const
2103 const size_t polygonCount = mMesh->polygonCount();
2104 const int subdivisionCount =
2105 polygonCount < SubTask::POLYGON_LIMIT ? evalSubdivisionCount(prim) : 0;
2107 if (subdivisionCount <= 0) {
2108 voxelizeTriangle(prim, data, mInterrupter);
2110 spawnTasks(prim, *mDataTable, subdivisionCount, polygonCount, mInterrupter);
2114 static void spawnTasks(
2115 const Triangle& mainPrim,
2116 DataTable& dataTable,
2117 int subdivisionCount,
2118 size_t polygonCount,
2119 Interrupter*
const interrupter)
2121 subdivisionCount -= 1;
2124 tbb::task_group tasks;
2126 const Vec3d ac = (mainPrim.a + mainPrim.c) * 0.5;
2127 const Vec3d bc = (mainPrim.b + mainPrim.c) * 0.5;
2128 const Vec3d ab = (mainPrim.a + mainPrim.b) * 0.5;
2131 prim.index = mainPrim.index;
2133 prim.a = mainPrim.a;
2136 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount, interrupter));
2141 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount, interrupter));
2144 prim.b = mainPrim.b;
2146 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount, interrupter));
2150 prim.c = mainPrim.c;
2151 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount, interrupter));
2156 static void voxelizeTriangle(
const Triangle& prim, VoxelizationDataType& data, Interrupter*
const interrupter)
2158 std::deque<Coord> coordList;
2162 coordList.push_back(ijk);
2167 updateDistance(ijk, prim, data);
2169 unsigned char primId = data.getNewPrimId();
2170 data.primIdAcc.setValueOnly(ijk, primId);
2172 while (!coordList.empty()) {
2173 if (interrupter && interrupter->wasInterrupted()) {
2174 thread::cancelGroupExecution();
2177 for (Int32 pass = 0; pass < 1048576 && !coordList.empty(); ++pass) {
2178 ijk = coordList.back();
2179 coordList.pop_back();
2181 for (Int32 i = 0; i < 26; ++i) {
2182 nijk = ijk + util::COORD_OFFSETS[i];
2183 if (primId != data.primIdAcc.getValue(nijk)) {
2184 data.primIdAcc.setValueOnly(nijk, primId);
2185 if(updateDistance(nijk, prim, data)) coordList.push_back(nijk);
2192 static bool updateDistance(
const Coord& ijk,
const Triangle& prim, VoxelizationDataType& data)
2194 Vec3d uvw, voxelCenter(ijk[0], ijk[1], ijk[2]);
2203 if (std::isnan(dist))
2206 const ValueType oldDist = data.distAcc.getValue(ijk);
2208 if (dist < oldDist) {
2209 data.distAcc.setValue(ijk, dist);
2210 data.indexAcc.setValue(ijk, prim.index);
2214 data.indexAcc.setValueOnly(ijk,
std::min(prim.index, data.indexAcc.getValue(ijk)));
2217 return !(dist > 0.75);
2220 DataTable *
const mDataTable;
2222 Interrupter *
const mInterrupter;
2229 template<
typename TreeType>
2230 struct DiffLeafNodeMask
2232 using AccessorType =
typename tree::ValueAccessor<TreeType>;
2233 using LeafNodeType =
typename TreeType::LeafNodeType;
2236 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
2238 DiffLeafNodeMask(
const TreeType& rhsTree,
2239 std::vector<BoolLeafNodeType*>& lhsNodes)
2240 : mRhsTree(&rhsTree), mLhsNodes(lhsNodes.empty() ? nullptr : &lhsNodes[0])
2244 void operator()(
const tbb::blocked_range<size_t>& range)
const {
2246 tree::ValueAccessor<const TreeType> acc(*mRhsTree);
2248 for (
size_t n = range.begin(), N = range.end(); n <
N; ++
n) {
2250 BoolLeafNodeType* lhsNode = mLhsNodes[
n];
2251 const LeafNodeType* rhsNode = acc.probeConstLeaf(lhsNode->origin());
2253 if (rhsNode) lhsNode->topologyDifference(*rhsNode,
false);
2258 TreeType
const *
const mRhsTree;
2259 BoolLeafNodeType **
const mLhsNodes;
2263 template<
typename LeafNodeTypeA,
typename LeafNodeTypeB>
2264 struct UnionValueMasks
2266 UnionValueMasks(std::vector<LeafNodeTypeA*>& nodesA, std::vector<LeafNodeTypeB*>& nodesB)
2267 : mNodesA(nodesA.empty() ? nullptr : &nodesA[0])
2268 , mNodesB(nodesB.empty() ? nullptr : &nodesB[0])
2272 void operator()(
const tbb::blocked_range<size_t>& range)
const {
2273 for (
size_t n = range.begin(), N = range.end(); n <
N; ++
n) {
2274 mNodesA[
n]->topologyUnion(*mNodesB[n]);
2279 LeafNodeTypeA **
const mNodesA;
2280 LeafNodeTypeB **
const mNodesB;
2284 template<
typename TreeType>
2285 struct ConstructVoxelMask
2287 using LeafNodeType =
typename TreeType::LeafNodeType;
2290 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
2292 ConstructVoxelMask(BoolTreeType& maskTree,
const TreeType& tree,
2293 std::vector<LeafNodeType*>& nodes)
2295 , mNodes(nodes.empty() ? nullptr : &nodes[0])
2296 , mLocalMaskTree(false)
2297 , mMaskTree(&maskTree)
2301 ConstructVoxelMask(ConstructVoxelMask& rhs,
tbb::split)
2303 , mNodes(rhs.mNodes)
2304 , mLocalMaskTree(false)
2305 , mMaskTree(&mLocalMaskTree)
2309 void operator()(
const tbb::blocked_range<size_t>& range)
2311 using Iterator =
typename LeafNodeType::ValueOnCIter;
2313 tree::ValueAccessor<const TreeType> acc(*mTree);
2314 tree::ValueAccessor<BoolTreeType> maskAcc(*mMaskTree);
2316 Coord ijk, nijk, localCorod;
2319 for (
size_t n = range.begin(); n != range.end(); ++
n) {
2321 LeafNodeType& node = *mNodes[
n];
2323 CoordBBox bbox = node.getNodeBoundingBox();
2326 BoolLeafNodeType& maskNode = *maskAcc.touchLeaf(node.origin());
2328 for (Iterator it = node.cbeginValueOn(); it; ++it) {
2329 ijk = it.getCoord();
2332 localCorod = LeafNodeType::offsetToLocalCoord(pos);
2334 if (localCorod[2] <
int(LeafNodeType::DIM - 1)) {
2336 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2338 nijk = ijk.offsetBy(0, 0, 1);
2339 if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2342 if (localCorod[2] > 0) {
2344 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2346 nijk = ijk.offsetBy(0, 0, -1);
2347 if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2350 if (localCorod[1] <
int(LeafNodeType::DIM - 1)) {
2351 npos = pos + LeafNodeType::DIM;
2352 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2354 nijk = ijk.offsetBy(0, 1, 0);
2355 if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2358 if (localCorod[1] > 0) {
2359 npos = pos - LeafNodeType::DIM;
2360 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2362 nijk = ijk.offsetBy(0, -1, 0);
2363 if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2366 if (localCorod[0] <
int(LeafNodeType::DIM - 1)) {
2367 npos = pos + LeafNodeType::DIM * LeafNodeType::DIM;
2368 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2370 nijk = ijk.offsetBy(1, 0, 0);
2371 if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2374 if (localCorod[0] > 0) {
2375 npos = pos - LeafNodeType::DIM * LeafNodeType::DIM;
2376 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2378 nijk = ijk.offsetBy(-1, 0, 0);
2379 if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2385 void join(ConstructVoxelMask& rhs) { mMaskTree->merge(*rhs.mMaskTree); }
2388 TreeType
const *
const mTree;
2389 LeafNodeType **
const mNodes;
2391 BoolTreeType mLocalMaskTree;
2392 BoolTreeType *
const mMaskTree;
2397 template<
typename TreeType,
typename MeshDataAdapter>
2398 struct ExpandNarrowband
2401 using LeafNodeType =
typename TreeType::LeafNodeType;
2402 using NodeMaskType =
typename LeafNodeType::NodeMaskType;
2404 using Int32LeafNodeType =
typename Int32TreeType::LeafNodeType;
2406 using BoolLeafNodeType =
typename BoolTreeType::LeafNodeType;
2413 Fragment() : idx(0),
x(0),
y(0), z(0), dist(0.0) {}
2415 Fragment(Int32 idx_, Int32 x_, Int32 y_, Int32 z_,
ValueType dist_)
2416 : idx(idx_),
x(x_),
y(y_), z(z_), dist(dist_)
2420 bool operator<(
const Fragment& rhs)
const {
return idx < rhs.idx; }
2426 std::vector<BoolLeafNodeType*>& maskNodes,
2427 BoolTreeType& maskTree,
2429 Int32TreeType& indexTree,
2434 : mMaskNodes(maskNodes.empty() ? nullptr : &maskNodes[0])
2435 , mMaskTree(&maskTree)
2436 , mDistTree(&distTree)
2437 , mIndexTree(&indexTree)
2439 , mNewMaskTree(false)
2441 , mUpdatedDistNodes()
2443 , mUpdatedIndexNodes()
2444 , mExteriorBandWidth(exteriorBandWidth)
2445 , mInteriorBandWidth(interiorBandWidth)
2446 , mVoxelSize(voxelSize)
2450 ExpandNarrowband(
const ExpandNarrowband& rhs,
tbb::split)
2451 : mMaskNodes(rhs.mMaskNodes)
2452 , mMaskTree(rhs.mMaskTree)
2453 , mDistTree(rhs.mDistTree)
2454 , mIndexTree(rhs.mIndexTree)
2456 , mNewMaskTree(false)
2458 , mUpdatedDistNodes()
2460 , mUpdatedIndexNodes()
2461 , mExteriorBandWidth(rhs.mExteriorBandWidth)
2462 , mInteriorBandWidth(rhs.mInteriorBandWidth)
2463 , mVoxelSize(rhs.mVoxelSize)
2467 void join(ExpandNarrowband& rhs)
2469 mDistNodes.insert(mDistNodes.end(), rhs.mDistNodes.begin(), rhs.mDistNodes.end());
2470 mIndexNodes.insert(mIndexNodes.end(), rhs.mIndexNodes.begin(), rhs.mIndexNodes.end());
2472 mUpdatedDistNodes.insert(mUpdatedDistNodes.end(),
2473 rhs.mUpdatedDistNodes.begin(), rhs.mUpdatedDistNodes.end());
2475 mUpdatedIndexNodes.insert(mUpdatedIndexNodes.end(),
2476 rhs.mUpdatedIndexNodes.begin(), rhs.mUpdatedIndexNodes.end());
2478 mNewMaskTree.merge(rhs.mNewMaskTree);
2481 void operator()(
const tbb::blocked_range<size_t>& range)
2483 tree::ValueAccessor<BoolTreeType> newMaskAcc(mNewMaskTree);
2484 tree::ValueAccessor<TreeType> distAcc(*mDistTree);
2485 tree::ValueAccessor<Int32TreeType> indexAcc(*mIndexTree);
2487 std::vector<Fragment> fragments;
2488 fragments.reserve(256);
2490 std::unique_ptr<LeafNodeType> newDistNodePt;
2491 std::unique_ptr<Int32LeafNodeType> newIndexNodePt;
2493 for (
size_t n = range.begin(), N = range.end(); n <
N; ++
n) {
2495 BoolLeafNodeType& maskNode = *mMaskNodes[
n];
2496 if (maskNode.isEmpty())
continue;
2500 const Coord& origin = maskNode.origin();
2502 LeafNodeType * distNodePt = distAcc.probeLeaf(origin);
2503 Int32LeafNodeType * indexNodePt = indexAcc.probeLeaf(origin);
2507 bool usingNewNodes =
false;
2509 if (!distNodePt && !indexNodePt) {
2511 const ValueType backgroundDist = distAcc.getValue(origin);
2513 if (!newDistNodePt.get() && !newIndexNodePt.get()) {
2514 newDistNodePt.reset(
new LeafNodeType(origin, backgroundDist));
2515 newIndexNodePt.reset(
new Int32LeafNodeType(origin, indexAcc.getValue(origin)));
2518 if ((backgroundDist <
ValueType(0.0)) !=
2519 (newDistNodePt->getValue(0) <
ValueType(0.0))) {
2520 newDistNodePt->buffer().fill(backgroundDist);
2523 newDistNodePt->setOrigin(origin);
2524 newIndexNodePt->setOrigin(origin);
2527 distNodePt = newDistNodePt.get();
2528 indexNodePt = newIndexNodePt.get();
2530 usingNewNodes =
true;
2537 for (
typename BoolLeafNodeType::ValueOnIter it = maskNode.beginValueOn(); it; ++it) {
2538 bbox.expand(it.getCoord());
2543 gatherFragments(fragments, bbox, distAcc, indexAcc);
2548 bbox = maskNode.getNodeBoundingBox();
2550 bool updatedLeafNodes =
false;
2552 for (
typename BoolLeafNodeType::ValueOnIter it = maskNode.beginValueOn(); it; ++it) {
2554 const Coord ijk = it.getCoord();
2556 if (updateVoxel(ijk, 5, fragments, *distNodePt, *indexNodePt, &updatedLeafNodes)) {
2558 for (Int32 i = 0; i < 6; ++i) {
2559 const Coord nijk = ijk + util::COORD_OFFSETS[i];
2560 if (bbox.isInside(nijk)) {
2561 mask.setOn(BoolLeafNodeType::coordToOffset(nijk));
2563 newMaskAcc.setValueOn(nijk);
2567 for (Int32 i = 6; i < 26; ++i) {
2568 const Coord nijk = ijk + util::COORD_OFFSETS[i];
2569 if (bbox.isInside(nijk)) {
2570 mask.setOn(BoolLeafNodeType::coordToOffset(nijk));
2576 if (updatedLeafNodes) {
2579 mask -= indexNodePt->getValueMask();
2581 for (
typename NodeMaskType::OnIterator it = mask.beginOn(); it; ++it) {
2583 const Index pos = it.pos();
2584 const Coord ijk = maskNode.origin() + LeafNodeType::offsetToLocalCoord(pos);
2586 if (updateVoxel(ijk, 6, fragments, *distNodePt, *indexNodePt)) {
2587 for (Int32 i = 0; i < 6; ++i) {
2588 newMaskAcc.setValueOn(ijk + util::COORD_OFFSETS[i]);
2594 if (usingNewNodes) {
2595 newDistNodePt->topologyUnion(*newIndexNodePt);
2596 mDistNodes.push_back(newDistNodePt.release());
2597 mIndexNodes.push_back(newIndexNodePt.release());
2599 mUpdatedDistNodes.push_back(distNodePt);
2600 mUpdatedIndexNodes.push_back(indexNodePt);
2608 BoolTreeType& newMaskTree() {
return mNewMaskTree; }
2610 std::vector<LeafNodeType*>& newDistNodes() {
return mDistNodes; }
2611 std::vector<LeafNodeType*>& updatedDistNodes() {
return mUpdatedDistNodes; }
2613 std::vector<Int32LeafNodeType*>& newIndexNodes() {
return mIndexNodes; }
2614 std::vector<Int32LeafNodeType*>& updatedIndexNodes() {
return mUpdatedIndexNodes; }
2620 gatherFragments(std::vector<Fragment>& fragments,
const CoordBBox& bbox,
2621 tree::ValueAccessor<TreeType>& distAcc, tree::ValueAccessor<Int32TreeType>& indexAcc)
2624 const Coord nodeMin = bbox.min() & ~(LeafNodeType::DIM - 1);
2625 const Coord nodeMax = bbox.max() & ~(LeafNodeType::DIM - 1);
2630 for (ijk[0] = nodeMin[0]; ijk[0] <= nodeMax[0]; ijk[0] += LeafNodeType::DIM) {
2631 for (ijk[1] = nodeMin[1]; ijk[1] <= nodeMax[1]; ijk[1] += LeafNodeType::DIM) {
2632 for (ijk[2] = nodeMin[2]; ijk[2] <= nodeMax[2]; ijk[2] += LeafNodeType::DIM) {
2633 if (LeafNodeType* distleaf = distAcc.probeLeaf(ijk)) {
2636 ijk.offsetBy(LeafNodeType::DIM - 1));
2637 gatherFragments(fragments, region, *distleaf, *indexAcc.probeLeaf(ijk));
2643 std::sort(fragments.begin(), fragments.end());
2647 gatherFragments(std::vector<Fragment>& fragments,
const CoordBBox& bbox,
2648 const LeafNodeType& distLeaf,
const Int32LeafNodeType& idxLeaf)
const
2650 const typename LeafNodeType::NodeMaskType& mask = distLeaf.getValueMask();
2651 const ValueType* distData = distLeaf.buffer().data();
2652 const Int32* idxData = idxLeaf.buffer().data();
2654 for (
int x = bbox.min()[0];
x <= bbox.max()[0]; ++
x) {
2655 const Index xPos = (
x & (LeafNodeType::DIM - 1u)) << (2 * LeafNodeType::LOG2DIM);
2656 for (
int y = bbox.min()[1];
y <= bbox.max()[1]; ++
y) {
2657 const Index yPos = xPos + ((
y & (LeafNodeType::DIM - 1u)) << LeafNodeType::LOG2DIM);
2658 for (
int z = bbox.min()[2]; z <= bbox.max()[2]; ++
z) {
2659 const Index pos = yPos + (z & (LeafNodeType::DIM - 1u));
2660 if (mask.isOn(pos)) {
2661 fragments.push_back(Fragment(idxData[pos],
x,
y,z,
std::abs(distData[pos])));
2671 computeDistance(
const Coord& ijk,
const Int32 manhattanLimit,
2672 const std::vector<Fragment>& fragments, Int32& closestPrimIdx)
const
2674 Vec3d a,
b,
c, uvw, voxelCenter(ijk[0], ijk[1], ijk[2]);
2678 for (
size_t n = 0, N = fragments.size(); n <
N; ++
n) {
2680 const Fragment& fragment = fragments[
n];
2681 if (lastIdx == fragment.idx)
continue;
2687 const Int32 manhattan = dx + dy + dz;
2688 if (manhattan > manhattanLimit)
continue;
2690 lastIdx = fragment.idx;
2692 const size_t polygon = size_t(lastIdx);
2694 mMesh->getIndexSpacePoint(polygon, 0, a);
2695 mMesh->getIndexSpacePoint(polygon, 1, b);
2696 mMesh->getIndexSpacePoint(polygon, 2, c);
2698 primDist = (voxelCenter -
2702 if (4 == mMesh->vertexCount(polygon)) {
2704 mMesh->getIndexSpacePoint(polygon, 3, b);
2707 a, b, c, voxelCenter, uvw)).lengthSqr();
2709 if (tmpDist < primDist) primDist = tmpDist;
2712 if (primDist < dist) {
2714 closestPrimIdx = lastIdx;
2724 updateVoxel(
const Coord& ijk,
const Int32 manhattanLimit,
2725 const std::vector<Fragment>& fragments,
2726 LeafNodeType& distLeaf, Int32LeafNodeType& idxLeaf,
bool* updatedLeafNodes =
nullptr)
2728 Int32 closestPrimIdx = 0;
2729 const ValueType distance = computeDistance(ijk, manhattanLimit, fragments, closestPrimIdx);
2731 const Index pos = LeafNodeType::coordToOffset(ijk);
2732 const bool inside = distLeaf.getValue(pos) <
ValueType(0.0);
2734 bool activateNeighbourVoxels =
false;
2736 if (!inside && distance < mExteriorBandWidth) {
2737 if (updatedLeafNodes) *updatedLeafNodes =
true;
2738 activateNeighbourVoxels = (distance + mVoxelSize) < mExteriorBandWidth;
2739 distLeaf.setValueOnly(pos, distance);
2740 idxLeaf.setValueOn(pos, closestPrimIdx);
2741 }
else if (inside && distance < mInteriorBandWidth) {
2742 if (updatedLeafNodes) *updatedLeafNodes =
true;
2743 activateNeighbourVoxels = (distance + mVoxelSize) < mInteriorBandWidth;
2744 distLeaf.setValueOnly(pos, -distance);
2745 idxLeaf.setValueOn(pos, closestPrimIdx);
2748 return activateNeighbourVoxels;
2753 BoolLeafNodeType **
const mMaskNodes;
2754 BoolTreeType *
const mMaskTree;
2755 TreeType *
const mDistTree;
2756 Int32TreeType *
const mIndexTree;
2760 BoolTreeType mNewMaskTree;
2762 std::vector<LeafNodeType*> mDistNodes, mUpdatedDistNodes;
2763 std::vector<Int32LeafNodeType*> mIndexNodes, mUpdatedIndexNodes;
2765 const ValueType mExteriorBandWidth, mInteriorBandWidth, mVoxelSize;
2769 template<
typename TreeType>
2771 using LeafNodeType =
typename TreeType::LeafNodeType;
2773 AddNodes(TreeType& tree, std::vector<LeafNodeType*>& nodes)
2774 : mTree(&tree) , mNodes(&nodes)
2778 void operator()()
const {
2779 tree::ValueAccessor<TreeType> acc(*mTree);
2780 std::vector<LeafNodeType*>& nodes = *mNodes;
2781 for (
size_t n = 0, N = nodes.size(); n <
N; ++
n) {
2782 acc.addLeaf(nodes[n]);
2786 TreeType *
const mTree;
2787 std::vector<LeafNodeType*> *
const mNodes;
2791 template<
typename TreeType,
typename Int32TreeType,
typename BoolTreeType,
typename MeshDataAdapter>
2795 Int32TreeType& indexTree,
2796 BoolTreeType& maskTree,
2797 std::vector<typename BoolTreeType::LeafNodeType*>& maskNodes,
2803 ExpandNarrowband<TreeType, MeshDataAdapter> expandOp(maskNodes, maskTree,
2804 distTree, indexTree, mesh, exteriorBandWidth, interiorBandWidth, voxelSize);
2806 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, maskNodes.size()), expandOp);
2808 tbb::parallel_for(tbb::blocked_range<size_t>(0, expandOp.updatedIndexNodes().size()),
2809 UnionValueMasks<typename TreeType::LeafNodeType, typename Int32TreeType::LeafNodeType>(
2810 expandOp.updatedDistNodes(), expandOp.updatedIndexNodes()));
2812 tbb::task_group tasks;
2813 tasks.run(AddNodes<TreeType>(distTree, expandOp.newDistNodes()));
2814 tasks.run(AddNodes<Int32TreeType>(indexTree, expandOp.newIndexNodes()));
2818 maskTree.merge(expandOp.newMaskTree());
2826 template<
typename TreeType>
2827 struct TransformValues
2829 using LeafNodeType =
typename TreeType::LeafNodeType;
2832 TransformValues(std::vector<LeafNodeType*>& nodes,
2835 , mVoxelSize(voxelSize)
2836 , mUnsigned(unsignedDist)
2840 void operator()(
const tbb::blocked_range<size_t>& range)
const {
2842 typename LeafNodeType::ValueOnIter iter;
2844 const bool udf = mUnsigned;
2845 const ValueType w[2] = { -mVoxelSize, mVoxelSize };
2847 for (
size_t n = range.begin(), N = range.end(); n <
N; ++
n) {
2849 for (iter = mNodes[n]->beginValueOn(); iter; ++iter) {
2857 LeafNodeType * *
const mNodes;
2859 const bool mUnsigned;
2864 template<
typename TreeType>
2865 struct InactivateValues
2867 using LeafNodeType =
typename TreeType::LeafNodeType;
2870 InactivateValues(std::vector<LeafNodeType*>& nodes,
2872 : mNodes(nodes.empty() ? nullptr : &nodes[0])
2873 , mExBandWidth(exBandWidth)
2874 , mInBandWidth(inBandWidth)
2878 void operator()(
const tbb::blocked_range<size_t>& range)
const {
2880 typename LeafNodeType::ValueOnIter iter;
2884 for (
size_t n = range.begin(), N = range.end(); n <
N; ++
n) {
2886 for (iter = mNodes[n]->beginValueOn(); iter; ++iter) {
2890 const bool inside = val <
ValueType(0.0);
2892 if (inside && !(val > inVal)) {
2895 }
else if (!inside && !(val < exVal)) {
2904 LeafNodeType * *
const mNodes;
2905 const ValueType mExBandWidth, mInBandWidth;
2909 template<
typename TreeType>
2912 using LeafNodeType =
typename TreeType::LeafNodeType;
2915 OffsetValues(std::vector<LeafNodeType*>& nodes,
ValueType offset)
2916 : mNodes(nodes.empty() ? nullptr : &nodes[0]),
mOffset(offset)
2920 void operator()(
const tbb::blocked_range<size_t>& range)
const {
2924 for (
size_t n = range.begin(), N = range.end(); n <
N; ++
n) {
2926 typename LeafNodeType::ValueOnIter iter = mNodes[
n]->beginValueOn();
2928 for (; iter; ++iter) {
2936 LeafNodeType * *
const mNodes;
2941 template<
typename TreeType>
2944 using LeafNodeType =
typename TreeType::LeafNodeType;
2947 Renormalize(
const TreeType& tree,
const std::vector<LeafNodeType*>& nodes,
2950 , mNodes(nodes.empty() ? nullptr : &nodes[0])
2952 , mVoxelSize(voxelSize)
2956 void operator()(
const tbb::blocked_range<size_t>& range)
const
2958 using Vec3Type = math::Vec3<ValueType>;
2960 tree::ValueAccessor<const TreeType> acc(*mTree);
2967 for (
size_t n = range.begin(), N = range.end(); n <
N; ++
n) {
2971 typename LeafNodeType::ValueOnCIter iter = mNodes[
n]->cbeginValueOn();
2972 for (; iter; ++iter) {
2976 ijk = iter.getCoord();
2978 up[0] = acc.getValue(ijk.offsetBy(1, 0, 0)) - phi0;
2979 up[1] = acc.getValue(ijk.offsetBy(0, 1, 0)) - phi0;
2980 up[2] = acc.getValue(ijk.offsetBy(0, 0, 1)) - phi0;
2982 down[0] = phi0 - acc.getValue(ijk.offsetBy(-1, 0, 0));
2983 down[1] = phi0 - acc.getValue(ijk.offsetBy(0, -1, 0));
2984 down[2] = phi0 - acc.getValue(ijk.offsetBy(0, 0, -1));
2991 bufferData[iter.pos()] = phi0 - dx * S * diff;
2997 TreeType
const *
const mTree;
2998 LeafNodeType
const *
const *
const mNodes;
3005 template<
typename TreeType>
3008 using LeafNodeType =
typename TreeType::LeafNodeType;
3011 MinCombine(std::vector<LeafNodeType*>& nodes,
const ValueType*
buffer)
3012 : mNodes(nodes.empty() ? nullptr : &nodes[0]), mBuffer(buffer)
3016 void operator()(
const tbb::blocked_range<size_t>& range)
const {
3018 for (
size_t n = range.begin(), N = range.end(); n <
N; ++
n) {
3022 typename LeafNodeType::ValueOnIter iter = mNodes[
n]->beginValueOn();
3024 for (; iter; ++iter) {
3026 val =
std::min(val, bufferData[iter.pos()]);
3032 LeafNodeType * *
const mNodes;
3047 template <
typename FloatTreeT>
3051 using ConnectivityTable = mesh_to_volume_internal::LeafNodeConnectivityTable<FloatTreeT>;
3056 ConnectivityTable nodeConnectivity(tree);
3058 std::vector<size_t> zStartNodes, yStartNodes, xStartNodes;
3063 for (
size_t n = 0; n < nodeConnectivity.size(); ++
n) {
3064 if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevX()[
n]) {
3065 xStartNodes.push_back(n);
3068 if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevY()[
n]) {
3069 yStartNodes.push_back(n);
3072 if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevZ()[
n]) {
3073 zStartNodes.push_back(n);
3077 using SweepingOp = mesh_to_volume_internal::SweepExteriorSign<FloatTreeT>;
3091 const size_t numLeafNodes = nodeConnectivity.size();
3094 std::unique_ptr<bool[]> changedNodeMaskA{
new bool[numLeafNodes]};
3095 std::unique_ptr<bool[]> changedNodeMaskB{
new bool[numLeafNodes]};
3096 std::unique_ptr<bool[]> changedVoxelMask{
new bool[numVoxels]};
3098 mesh_to_volume_internal::fillArray(changedNodeMaskA.get(),
true, numLeafNodes);
3099 mesh_to_volume_internal::fillArray(changedNodeMaskB.get(),
false, numLeafNodes);
3100 mesh_to_volume_internal::fillArray(changedVoxelMask.get(),
false, numVoxels);
3102 const tbb::blocked_range<size_t> nodeRange(0, numLeafNodes);
3104 bool nodesUpdated =
false;
3109 tbb::parallel_for(nodeRange, mesh_to_volume_internal::SeedFillExteriorSign<FloatTreeT>(
3110 nodeConnectivity.nodes(), changedNodeMaskA.get()));
3119 nodeConnectivity, changedNodeMaskA.get(), changedNodeMaskB.get(),
3120 changedVoxelMask.get()));
3124 changedNodeMaskA.swap(changedNodeMaskB);
3126 nodesUpdated =
false;
3127 for (
size_t n = 0; n < numLeafNodes; ++
n) {
3128 nodesUpdated |= changedNodeMaskA[
n];
3129 if (nodesUpdated)
break;
3135 tbb::parallel_for(nodeRange, mesh_to_volume_internal::SyncVoxelMask<FloatTreeT>(
3136 nodeConnectivity.nodes(), changedNodeMaskA.get(), changedVoxelMask.get()));
3138 }
while (nodesUpdated);
3145 template <
typename T, Index Log2Dim,
typename InteriorTest>
3168 const auto DIM = leafNode.
DIM;
3171 std::vector<VoxelState> voxelState(
SIZE, NOT_VISITED);
3173 std::vector<std::pair<Index, VoxelState>> offsetStack;
3174 offsetStack.reserve(
SIZE);
3176 for (
Index offset=0; offset<
SIZE; offset++) {
3183 }
else if (voxelState[offset] == NOT_VISITED) {
3187 if (interiorTest(coord)){
3190 offsetStack.push_back({
offset, POSITIVE});
3191 voxelState[
offset] = POSITIVE;
3193 offsetStack.push_back({
offset, NEGATIVE});
3194 voxelState[
offset] = NEGATIVE;
3197 while(!offsetStack.empty()){
3199 auto [off,
state] = offsetStack[offsetStack.size()-1];
3200 offsetStack.pop_back();
3202 if (
state == NEGATIVE) {
3209 for (
int dim=2; dim>=0; dim--){
3210 for (
int i = -1; i <=1; ++(++i)){
3211 int dimIdx = (off >> dim * Log2Dim) % DIM;
3212 auto neighOff = off + (1 << dim * Log2Dim) * i;
3214 (dimIdx < (
int)DIM - 1) &&
3215 (voxelState[neighOff] == NOT_VISITED)) {
3220 offsetStack.push_back({neighOff,
state});
3221 voxelState[neighOff] =
state;
3251 template <
typename FloatTreeT,
typename InteriorTest>
3256 "InteriorTest has to be a function `Coord -> bool`!");
3257 static_assert(std::is_copy_constructible_v<InteriorTest>,
3258 "InteriorTest has to be copyable!");
3260 using LeafT =
typename FloatTreeT::LeafNodeType;
3264 auto op = [interiorTest](
auto& node) {
3265 using Node = std::decay_t<decltype(node)>;
3267 if constexpr (std::is_same_v<Node, LeafT>) {
3269 for (
auto iter = node.beginValueAll(); iter; ++iter) {
3270 if (!interiorTest(iter.getCoord())) {
3271 iter.setValue(-*iter);
3276 for (
auto iter = node.beginChildOff(); iter; ++iter) {
3277 if (!interiorTest(iter.getCoord())) {
3278 iter.setValue(-*iter);
3284 openvdb::tree::NodeManager nodes(tree);
3285 nodes.foreachBottomUp(
op);
3290 auto op = [interiorTest](
auto& node) {
3291 using Node = std::decay_t<decltype(node)>;
3293 if constexpr (std::is_same_v<Node, LeafT>) {
3295 LeafT& leaf =
static_cast<LeafT&
>(node);
3300 for (
auto iter = node.beginChildOff(); iter; ++iter) {
3301 if (!interiorTest(iter.getCoord())) {
3302 iter.setValue(-*iter);
3308 openvdb::tree::NodeManager nodes(tree);
3309 nodes.foreachBottomUp(
op);
3316 template <
typename Gr
idType,
typename MeshDataAdapter,
typename Interrupter,
typename InteriorTest>
3317 typename GridType::Ptr
3319 Interrupter& interrupter,
3322 float exteriorBandWidth,
3323 float interiorBandWidth,
3326 InteriorTest interiorTest,
3329 using GridTypePtr =
typename GridType::Ptr;
3330 using TreeType =
typename GridType::TreeType;
3331 using LeafNodeType =
typename TreeType::LeafNodeType;
3335 using Int32TreeType =
typename Int32GridType::TreeType;
3344 distGrid->setTransform(transform.
copy());
3351 if (!std::isfinite(exteriorWidth) || std::isnan(interiorWidth)) {
3352 std::stringstream msg;
3353 msg <<
"Illegal narrow band width: exterior = " << exteriorWidth
3354 <<
", interior = " << interiorWidth;
3361 if (!std::isfinite(voxelSize) ||
math::isZero(voxelSize)) {
3362 std::stringstream msg;
3363 msg <<
"Illegal transform, voxel size = " << voxelSize;
3369 exteriorWidth *= voxelSize;
3373 interiorWidth *= voxelSize;
3381 Int32GridType* indexGrid =
nullptr;
3383 typename Int32GridType::Ptr temporaryIndexGrid;
3385 if (polygonIndexGrid) {
3386 indexGrid = polygonIndexGrid;
3389 indexGrid = temporaryIndexGrid.get();
3392 indexGrid->newTree();
3393 indexGrid->setTransform(transform.
copy());
3395 if (computeSignedDistanceField) {
3402 TreeType& distTree = distGrid->tree();
3403 Int32TreeType& indexTree = indexGrid->tree();
3411 using VoxelizationDataType = mesh_to_volume_internal::VoxelizationData<TreeType>;
3412 using DataTable = tbb::enumerable_thread_specific<typename VoxelizationDataType::Ptr>;
3416 mesh_to_volume_internal::VoxelizePolygons<TreeType, MeshDataAdapter, Interrupter>;
3418 const tbb::blocked_range<size_t> polygonRange(0, mesh.polygonCount());
3422 for (
typename DataTable::iterator i = data.begin(); i != data.end(); ++i) {
3423 VoxelizationDataType& dataItem = **i;
3424 mesh_to_volume_internal::combineData(
3425 distTree, indexTree, dataItem.distTree, dataItem.indexTree);
3431 if (interrupter.wasInterrupted(30))
return distGrid;
3438 if (computeSignedDistanceField) {
3441 if constexpr (std::is_same_v<InteriorTest, std::nullptr_t>) {
3443 (
void) interiorTest;
3450 bool signInitializedForEveryVoxel =
3452 !std::is_same_v<InteriorTest, std::nullptr_t> &&
3456 if (!signInitializedForEveryVoxel) {
3458 std::vector<LeafNodeType*> nodes;
3459 nodes.reserve(distTree.leafCount());
3460 distTree.getNodes(nodes);
3462 const tbb::blocked_range<size_t> nodeRange(0, nodes.size());
3465 mesh_to_volume_internal::ComputeIntersectingVoxelSign<TreeType, MeshDataAdapter>;
3469 if (interrupter.wasInterrupted(45))
return distGrid;
3472 if (removeIntersectingVoxels) {
3475 mesh_to_volume_internal::ValidateIntersectingVoxels<TreeType>(distTree, nodes));
3478 mesh_to_volume_internal::RemoveSelfIntersectingSurface<TreeType>(
3479 nodes, distTree, indexTree));
3487 if (interrupter.wasInterrupted(50))
return distGrid;
3489 if (distTree.activeVoxelCount() == 0) {
3491 distTree.root().setBackground(exteriorWidth,
false);
3497 std::vector<LeafNodeType*> nodes;
3498 nodes.reserve(distTree.leafCount());
3499 distTree.getNodes(nodes);
3502 mesh_to_volume_internal::TransformValues<TreeType>(
3503 nodes, voxelSize, !computeSignedDistanceField));
3507 if (computeSignedDistanceField) {
3508 distTree.root().setBackground(exteriorWidth,
false);
3514 if (interrupter.wasInterrupted(54))
return distGrid;
3523 if (interiorWidth > minBandWidth || exteriorWidth > minBandWidth) {
3526 BoolTreeType maskTree(
false);
3529 std::vector<LeafNodeType*> nodes;
3530 nodes.reserve(distTree.leafCount());
3531 distTree.getNodes(nodes);
3533 mesh_to_volume_internal::ConstructVoxelMask<TreeType>
op(maskTree, distTree, nodes);
3534 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), op);
3540 float progress = 54.0f, step = 0.0f;
3542 2.0 *
std::ceil((
std::max(interiorWidth, exteriorWidth) - minBandWidth) / voxelSize);
3544 if (estimated <
double(maxIterations)) {
3545 maxIterations = unsigned(estimated);
3546 step = 40.0f /
float(maxIterations);
3549 std::vector<typename BoolTreeType::LeafNodeType*> maskNodes;
3554 if (interrupter.wasInterrupted(
int(progress)))
return distGrid;
3556 const size_t maskNodeCount = maskTree.leafCount();
3557 if (maskNodeCount == 0)
break;
3560 maskNodes.reserve(maskNodeCount);
3561 maskTree.getNodes(maskNodes);
3563 const tbb::blocked_range<size_t>
range(0, maskNodes.size());
3566 mesh_to_volume_internal::DiffLeafNodeMask<TreeType>(distTree, maskNodes));
3568 mesh_to_volume_internal::expandNarrowband(distTree, indexTree, maskTree, maskNodes,
3569 mesh, exteriorWidth, interiorWidth, voxelSize);
3571 if ((++count) >= maxIterations)
break;
3576 if (interrupter.wasInterrupted(94))
return distGrid;
3578 if (!polygonIndexGrid) indexGrid->clear();
3586 if (computeSignedDistanceField && renormalizeValues) {
3588 std::vector<LeafNodeType*> nodes;
3589 nodes.reserve(distTree.leafCount());
3590 distTree.getNodes(nodes);
3592 std::unique_ptr<ValueType[]>
buffer{
new ValueType[LeafNodeType::SIZE * nodes.size()]};
3594 const ValueType offset =
ValueType(0.8 * voxelSize);
3597 mesh_to_volume_internal::OffsetValues<TreeType>(nodes, -offset));
3600 mesh_to_volume_internal::Renormalize<TreeType>(
3601 distTree, nodes,
buffer.get(), voxelSize));
3604 mesh_to_volume_internal::MinCombine<TreeType>(nodes,
buffer.get()));
3607 mesh_to_volume_internal::OffsetValues<TreeType>(
3608 nodes, offset - mesh_to_volume_internal::Tolerance<ValueType>::epsilon()));
3611 if (interrupter.wasInterrupted(99))
return distGrid;
3618 if (trimNarrowBand &&
std::min(interiorWidth, exteriorWidth) < voxelSize *
ValueType(4.0)) {
3620 std::vector<LeafNodeType*> nodes;
3621 nodes.reserve(distTree.leafCount());
3622 distTree.getNodes(nodes);
3625 mesh_to_volume_internal::InactivateValues<TreeType>(
3626 nodes, exteriorWidth, computeSignedDistanceField ? interiorWidth : exteriorWidth));
3629 distTree, exteriorWidth, computeSignedDistanceField ? -interiorWidth : -exteriorWidth);
3636 template <
typename Gr
idType,
typename MeshDataAdapter,
typename InteriorTest>
3637 typename GridType::Ptr
3641 float exteriorBandWidth,
3642 float interiorBandWidth,
3649 return meshToVolume<GridType>(nullInterrupter, mesh,
transform,
3650 exteriorBandWidth, interiorBandWidth,
flags, polygonIndexGrid);
3661 template<
typename Gr
idType,
typename Interrupter>
3663 typename GridType::Ptr>
::type
3665 Interrupter& interrupter,
3666 const openvdb::math::Transform& xform,
3667 const std::vector<Vec3s>& points,
3668 const std::vector<Vec3I>& triangles,
3669 const std::vector<Vec4I>& quads,
3672 bool unsignedDistanceField =
false)
3674 if (points.empty()) {
3678 const size_t numPoints = points.size();
3679 std::unique_ptr<Vec3s[]> indexSpacePoints{
new Vec3s[numPoints]};
3683 mesh_to_volume_internal::TransformPoints<Vec3s>(
3684 &points[0], indexSpacePoints.get(), xform));
3688 if (quads.empty()) {
3690 QuadAndTriangleDataAdapter<Vec3s, Vec3I>
3691 mesh(indexSpacePoints.get(), numPoints, &triangles[0], triangles.size());
3693 return meshToVolume<GridType>(
3694 interrupter, mesh, xform, exBandWidth, inBandWidth, conversionFlags);
3696 }
else if (triangles.empty()) {
3698 QuadAndTriangleDataAdapter<Vec3s, Vec4I>
3699 mesh(indexSpacePoints.get(), numPoints, &quads[0], quads.size());
3701 return meshToVolume<GridType>(
3702 interrupter, mesh, xform, exBandWidth, inBandWidth, conversionFlags);
3707 const size_t numPrimitives = triangles.size() + quads.size();
3708 std::unique_ptr<Vec4I[]> prims{
new Vec4I[numPrimitives]};
3710 for (
size_t n = 0, N = triangles.size(); n <
N; ++
n) {
3711 const Vec3I& triangle = triangles[
n];
3713 prim[0] = triangle[0];
3714 prim[1] = triangle[1];
3715 prim[2] = triangle[2];
3719 const size_t offset = triangles.size();
3720 for (
size_t n = 0, N = quads.size(); n <
N; ++
n) {
3721 prims[offset +
n] = quads[
n];
3724 QuadAndTriangleDataAdapter<Vec3s, Vec4I>
3725 mesh(indexSpacePoints.get(), numPoints, prims.get(), numPrimitives);
3727 return meshToVolume<GridType>(interrupter, mesh, xform,
3728 exBandWidth, inBandWidth, conversionFlags);
3734 template<
typename Gr
idType,
typename Interrupter>
3736 typename GridType::Ptr>
::type
3739 const math::Transform& ,
3740 const std::vector<Vec3s>& ,
3741 const std::vector<Vec3I>& ,
3742 const std::vector<Vec4I>& ,
3748 "mesh to volume conversion is supported only for scalar floating-point grids");
3758 template<
typename Gr
idType>
3759 typename GridType::Ptr
3761 const openvdb::math::Transform& xform,
3762 const std::vector<Vec3s>& points,
3763 const std::vector<Vec3I>& triangles,
3767 return meshToLevelSet<GridType>(nullInterrupter, xform,
points, triangles, halfWidth);
3771 template<
typename Gr
idType,
typename Interrupter>
3772 typename GridType::Ptr
3774 Interrupter& interrupter,
3775 const openvdb::math::Transform& xform,
3776 const std::vector<Vec3s>& points,
3777 const std::vector<Vec3I>& triangles,
3780 std::vector<Vec4I>
quads(0);
3781 return doMeshConversion<GridType>(interrupter, xform,
points, triangles,
quads,
3782 halfWidth, halfWidth);
3786 template<
typename Gr
idType>
3787 typename GridType::Ptr
3789 const openvdb::math::Transform& xform,
3790 const std::vector<Vec3s>& points,
3791 const std::vector<Vec4I>& quads,
3795 return meshToLevelSet<GridType>(nullInterrupter, xform,
points,
quads, halfWidth);
3799 template<
typename Gr
idType,
typename Interrupter>
3800 typename GridType::Ptr
3802 Interrupter& interrupter,
3803 const openvdb::math::Transform& xform,
3804 const std::vector<Vec3s>& points,
3805 const std::vector<Vec4I>& quads,
3808 std::vector<Vec3I> triangles(0);
3809 return doMeshConversion<GridType>(interrupter, xform,
points, triangles,
quads,
3810 halfWidth, halfWidth);
3814 template<
typename Gr
idType>
3815 typename GridType::Ptr
3817 const openvdb::math::Transform& xform,
3818 const std::vector<Vec3s>& points,
3819 const std::vector<Vec3I>& triangles,
3820 const std::vector<Vec4I>& quads,
3824 return meshToLevelSet<GridType>(
3825 nullInterrupter, xform,
points, triangles,
quads, halfWidth);
3829 template<
typename Gr
idType,
typename Interrupter>
3830 typename GridType::Ptr
3832 Interrupter& interrupter,
3833 const openvdb::math::Transform& xform,
3834 const std::vector<Vec3s>& points,
3835 const std::vector<Vec3I>& triangles,
3836 const std::vector<Vec4I>& quads,
3839 return doMeshConversion<GridType>(interrupter, xform,
points, triangles,
quads,
3840 halfWidth, halfWidth);
3844 template<
typename Gr
idType>
3845 typename GridType::Ptr
3847 const openvdb::math::Transform& xform,
3848 const std::vector<Vec3s>& points,
3849 const std::vector<Vec3I>& triangles,
3850 const std::vector<Vec4I>& quads,
3855 return meshToSignedDistanceField<GridType>(
3856 nullInterrupter, xform,
points, triangles,
quads, exBandWidth, inBandWidth);
3860 template<
typename Gr
idType,
typename Interrupter>
3861 typename GridType::Ptr
3863 Interrupter& interrupter,
3864 const openvdb::math::Transform& xform,
3865 const std::vector<Vec3s>& points,
3866 const std::vector<Vec3I>& triangles,
3867 const std::vector<Vec4I>& quads,
3871 return doMeshConversion<GridType>(interrupter, xform,
points, triangles,
3872 quads, exBandWidth, inBandWidth);
3876 template<
typename Gr
idType>
3877 typename GridType::Ptr
3879 const openvdb::math::Transform& xform,
3880 const std::vector<Vec3s>& points,
3881 const std::vector<Vec3I>& triangles,
3882 const std::vector<Vec4I>& quads,
3886 return meshToUnsignedDistanceField<GridType>(
3887 nullInterrupter, xform,
points, triangles,
quads, bandWidth);
3891 template<
typename Gr
idType,
typename Interrupter>
3892 typename GridType::Ptr
3894 Interrupter& interrupter,
3895 const openvdb::math::Transform& xform,
3896 const std::vector<Vec3s>& points,
3897 const std::vector<Vec3I>& triangles,
3898 const std::vector<Vec4I>& quads,
3901 return doMeshConversion<GridType>(interrupter, xform,
points, triangles,
quads,
3902 bandWidth, bandWidth,
true);
3910 inline std::ostream&
3913 ostr <<
"{[ " << rhs.
mXPrim <<
", " << rhs.
mXDist <<
"]";
3914 ostr <<
" [ " << rhs.
mYPrim <<
", " << rhs.
mYDist <<
"]";
3915 ostr <<
" [ " << rhs.
mZPrim <<
", " << rhs.
mZDist <<
"]}";
3920 inline MeshToVoxelEdgeData::EdgeData
3935 const std::vector<Vec3s>& pointList,
3936 const std::vector<Vec4I>& polygonList);
3941 inline void operator() (
const tbb::blocked_range<size_t> &range);
3951 template<
bool IsQuad>
3952 inline void voxelize(
const Primitive&);
3954 template<
bool IsQuad>
3955 inline bool evalPrimitive(
const Coord&,
const Primitive&);
3957 inline bool rayTriangleIntersection(
const Vec3d& origin,
const Vec3d& dir,
3958 const Vec3d& a,
const Vec3d& b,
const Vec3d& c,
double&
t);
3964 const std::vector<Vec3s>& mPointList;
3965 const std::vector<Vec4I>& mPolygonList;
3969 IntTreeT mLastPrimTree;
3976 const std::vector<Vec3s>& pointList,
3977 const std::vector<Vec4I>& polygonList)
3980 , mPointList(pointList)
3981 , mPolygonList(polygonList)
3983 , mLastPrimAccessor(mLastPrimTree)
3992 , mPointList(rhs.mPointList)
3993 , mPolygonList(rhs.mPolygonList)
3995 , mLastPrimAccessor(mLastPrimTree)
4004 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPolygonList.size()), *
this);
4006 (*this)(tbb::blocked_range<size_t>(0, mPolygonList.size()));
4015 using NodeChainType = RootNodeType::NodeChainType;
4016 static_assert(NodeChainType::Size > 1,
"expected tree height > 1");
4017 using InternalNodeType =
typename NodeChainType::template Get<1>;
4025 for ( ; leafIt; ++leafIt) {
4026 ijk = leafIt->origin();
4032 mAccessor.addLeaf(rhs.mAccessor.
probeLeaf(ijk));
4033 InternalNodeType* node = rhs.mAccessor.
getNode<InternalNodeType>();
4035 rhs.mAccessor.
clear();
4045 if (!lhsLeafPt->isValueOn(offset)) {
4046 lhsLeafPt->setValueOn(offset, rhsValue);
4078 for (
size_t n = range.begin(); n < range.end(); ++
n) {
4080 const Vec4I& verts = mPolygonList[
n];
4082 prim.index =
Int32(n);
4083 prim.a = Vec3d(mPointList[verts[0]]);
4084 prim.b = Vec3d(mPointList[verts[1]]);
4085 prim.c = Vec3d(mPointList[verts[2]]);
4088 prim.d = Vec3d(mPointList[verts[3]]);
4089 voxelize<true>(prim);
4091 voxelize<false>(prim);
4097 template<
bool IsQuad>
4099 MeshToVoxelEdgeData::GenEdgeData::voxelize(
const Primitive& prim)
4101 std::deque<Coord> coordList;
4105 coordList.push_back(ijk);
4107 evalPrimitive<IsQuad>(ijk, prim);
4109 while (!coordList.empty()) {
4111 ijk = coordList.back();
4112 coordList.pop_back();
4114 for (
Int32 i = 0; i < 26; ++i) {
4115 nijk = ijk + util::COORD_OFFSETS[i];
4117 if (prim.index != mLastPrimAccessor.getValue(nijk)) {
4118 mLastPrimAccessor.setValue(nijk, prim.index);
4119 if(evalPrimitive<IsQuad>(nijk, prim)) coordList.push_back(nijk);
4126 template<
bool IsQuad>
4128 MeshToVoxelEdgeData::GenEdgeData::evalPrimitive(
const Coord& ijk,
const Primitive& prim)
4130 Vec3d uvw, org(ijk[0], ijk[1], ijk[2]);
4131 bool intersecting =
false;
4135 mAccessor.probeValue(ijk, edgeData);
4138 double dist = (org -
4141 if (rayTriangleIntersection(org, Vec3d(1.0, 0.0, 0.0), prim.a, prim.c, prim.b, t)) {
4142 if (t < edgeData.mXDist) {
4143 edgeData.mXDist =
float(t);
4144 edgeData.mXPrim = prim.index;
4145 intersecting =
true;
4149 if (rayTriangleIntersection(org,
Vec3d(0.0, 1.0, 0.0), prim.a, prim.c, prim.b, t)) {
4150 if (t < edgeData.mYDist) {
4151 edgeData.mYDist =
float(t);
4152 edgeData.mYPrim = prim.index;
4153 intersecting =
true;
4157 if (rayTriangleIntersection(org,
Vec3d(0.0, 0.0, 1.0), prim.a, prim.c, prim.b, t)) {
4158 if (t < edgeData.mZDist) {
4159 edgeData.mZDist =
float(t);
4160 edgeData.mZPrim = prim.index;
4161 intersecting =
true;
4167 double secondDist = (org -
4170 if (secondDist < dist) dist = secondDist;
4172 if (rayTriangleIntersection(org,
Vec3d(1.0, 0.0, 0.0), prim.a, prim.d, prim.c, t)) {
4173 if (t < edgeData.mXDist) {
4174 edgeData.mXDist =
float(t);
4175 edgeData.mXPrim = prim.index;
4176 intersecting =
true;
4180 if (rayTriangleIntersection(org,
Vec3d(0.0, 1.0, 0.0), prim.a, prim.d, prim.c, t)) {
4181 if (t < edgeData.mYDist) {
4182 edgeData.mYDist =
float(t);
4183 edgeData.mYPrim = prim.index;
4184 intersecting =
true;
4188 if (rayTriangleIntersection(org,
Vec3d(0.0, 0.0, 1.0), prim.a, prim.d, prim.c, t)) {
4189 if (t < edgeData.mZDist) {
4190 edgeData.mZDist =
float(t);
4191 edgeData.mZPrim = prim.index;
4192 intersecting =
true;
4197 if (intersecting) mAccessor.setValue(ijk, edgeData);
4199 return (dist < 0.86602540378443861);
4204 MeshToVoxelEdgeData::GenEdgeData::rayTriangleIntersection(
4205 const Vec3d& origin,
const Vec3d& dir,
4206 const Vec3d& a,
const Vec3d& b,
const Vec3d& c,
4216 if (!(
std::abs(divisor) > 0.0))
return false;
4220 double inv_divisor = 1.0 /
divisor;
4222 double b1 = d.dot(s1) * inv_divisor;
4224 if (b1 < 0.0 || b1 > 1.0)
return false;
4227 double b2 = dir.
dot(s2) * inv_divisor;
4229 if (b2 < 0.0 || (b1 + b2) > 1.0)
return false;
4233 t = e2.dot(s2) * inv_divisor;
4234 return (t < 0.0) ?
false :
true;
4250 const std::vector<Vec3s>& pointList,
4251 const std::vector<Vec4I>& polygonList)
4265 std::vector<Vec3d>& points,
4266 std::vector<Index32>& primitives)
4276 point[0] = double(coord[0]) + data.
mXDist;
4277 point[1] = double(coord[1]);
4278 point[2] = double(coord[2]);
4280 points.push_back(point);
4281 primitives.push_back(data.
mXPrim);
4285 point[0] = double(coord[0]);
4286 point[1] = double(coord[1]) + data.
mYDist;
4287 point[2] = double(coord[2]);
4289 points.push_back(point);
4290 primitives.push_back(data.
mYPrim);
4294 point[0] = double(coord[0]);
4295 point[1] = double(coord[1]);
4296 point[2] = double(coord[2]) + data.
mZDist;
4298 points.push_back(point);
4299 primitives.push_back(data.
mZPrim);
4309 point[0] = double(coord[0]);
4310 point[1] = double(coord[1]) + data.
mYDist;
4311 point[2] = double(coord[2]);
4313 points.push_back(point);
4314 primitives.push_back(data.
mYPrim);
4318 point[0] = double(coord[0]);
4319 point[1] = double(coord[1]);
4320 point[2] = double(coord[2]) + data.
mZDist;
4322 points.push_back(point);
4323 primitives.push_back(data.
mZPrim);
4331 point[0] = double(coord[0]);
4332 point[1] = double(coord[1]) + data.
mYDist;
4333 point[2] = double(coord[2]);
4335 points.push_back(point);
4336 primitives.push_back(data.
mYPrim);
4345 point[0] = double(coord[0]) + data.
mXDist;
4346 point[1] = double(coord[1]);
4347 point[2] = double(coord[2]);
4349 points.push_back(point);
4350 primitives.push_back(data.
mXPrim);
4354 point[0] = double(coord[0]);
4355 point[1] = double(coord[1]) + data.
mYDist;
4356 point[2] = double(coord[2]);
4358 points.push_back(point);
4359 primitives.push_back(data.
mYPrim);
4369 point[0] = double(coord[0]) + data.
mXDist;
4370 point[1] = double(coord[1]);
4371 point[2] = double(coord[2]);
4373 points.push_back(point);
4374 primitives.push_back(data.
mXPrim);
4383 point[0] = double(coord[0]) + data.
mXDist;
4384 point[1] = double(coord[1]);
4385 point[2] = double(coord[2]);
4387 points.push_back(point);
4388 primitives.push_back(data.
mXPrim);
4392 point[0] = double(coord[0]);
4393 point[1] = double(coord[1]);
4394 point[2] = double(coord[2]) + data.
mZDist;
4396 points.push_back(point);
4397 primitives.push_back(data.
mZPrim);
4406 point[0] = double(coord[0]);
4407 point[1] = double(coord[1]);
4408 point[2] = double(coord[2]) + data.
mZDist;
4410 points.push_back(point);
4411 primitives.push_back(data.
mZPrim);
4417 template<
typename Gr
idType,
typename VecType>
4418 typename GridType::Ptr
4420 const openvdb::math::Transform& xform,
4427 points[0] =
Vec3s(pmin[0], pmin[1], pmin[2]);
4428 points[1] =
Vec3s(pmin[0], pmin[1], pmax[2]);
4429 points[2] =
Vec3s(pmax[0], pmin[1], pmax[2]);
4430 points[3] =
Vec3s(pmax[0], pmin[1], pmin[2]);
4431 points[4] =
Vec3s(pmin[0], pmax[1], pmin[2]);
4432 points[5] =
Vec3s(pmin[0], pmax[1], pmax[2]);
4433 points[6] =
Vec3s(pmax[0], pmax[1], pmax[2]);
4434 points[7] =
Vec3s(pmax[0], pmax[1], pmin[2]);
4437 faces[0] =
Vec4I(0, 1, 2, 3);
4438 faces[1] =
Vec4I(7, 6, 5, 4);
4439 faces[2] =
Vec4I(4, 5, 1, 0);
4440 faces[3] =
Vec4I(6, 7, 3, 2);
4441 faces[4] =
Vec4I(0, 3, 7, 4);
4442 faces[5] =
Vec4I(1, 5, 6, 2);
4446 return meshToVolume<GridType>(mesh, xform,
static_cast<float>(halfWidth), static_cast<float>(halfWidth));
4455 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
4457 #ifdef OPENVDB_INSTANTIATE_MESHTOVOLUME
4461 #define _FUNCTION(TreeT) \
4462 Grid<TreeT>::Ptr meshToVolume<Grid<TreeT>>(util::NullInterrupter&, \
4463 const QuadAndTriangleDataAdapter<Vec3s, Vec3I>&, const openvdb::math::Transform&, \
4464 float, float, int, Grid<TreeT>::ValueConverter<Int32>::Type*, std::nullptr_t, InteriorTestStrategy)
4468 #define _FUNCTION(TreeT) \
4469 Grid<TreeT>::Ptr meshToVolume<Grid<TreeT>>(util::NullInterrupter&, \
4470 const QuadAndTriangleDataAdapter<Vec3s, Vec4I>&, const openvdb::math::Transform&, \
4471 float, float, int, Grid<TreeT>::ValueConverter<Int32>::Type*, std::nullptr_t, InteriorTestStrategy)
4475 #define _FUNCTION(TreeT) \
4476 Grid<TreeT>::Ptr meshToLevelSet<Grid<TreeT>>(util::NullInterrupter&, \
4477 const openvdb::math::Transform&, const std::vector<Vec3s>&, const std::vector<Vec3I>&, \
4482 #define _FUNCTION(TreeT) \
4483 Grid<TreeT>::Ptr meshToLevelSet<Grid<TreeT>>(util::NullInterrupter&, \
4484 const openvdb::math::Transform&, const std::vector<Vec3s>&, const std::vector<Vec4I>&, \
4489 #define _FUNCTION(TreeT) \
4490 Grid<TreeT>::Ptr meshToLevelSet<Grid<TreeT>>(util::NullInterrupter&, \
4491 const openvdb::math::Transform&, const std::vector<Vec3s>&, \
4492 const std::vector<Vec3I>&, const std::vector<Vec4I>&, float)
4496 #define _FUNCTION(TreeT) \
4497 Grid<TreeT>::Ptr meshToSignedDistanceField<Grid<TreeT>>(util::NullInterrupter&, \
4498 const openvdb::math::Transform&, const std::vector<Vec3s>&, \
4499 const std::vector<Vec3I>&, const std::vector<Vec4I>&, float, float)
4503 #define _FUNCTION(TreeT) \
4504 Grid<TreeT>::Ptr meshToUnsignedDistanceField<Grid<TreeT>>(util::NullInterrupter&, \
4505 const openvdb::math::Transform&, const std::vector<Vec3s>&, \
4506 const std::vector<Vec3I>&, const std::vector<Vec4I>&, float)
4510 #define _FUNCTION(TreeT) \
4511 Grid<TreeT>::Ptr createLevelSetBox<Grid<TreeT>>(const math::BBox<Vec3s>&, \
4512 const openvdb::math::Transform&, float)
4516 #define _FUNCTION(TreeT) \
4517 Grid<TreeT>::Ptr createLevelSetBox<Grid<TreeT>>(const math::BBox<Vec3d>&, \
4518 const openvdb::math::Transform&, double)
4522 #define _FUNCTION(TreeT) \
4523 void traceExteriorBoundaries(TreeT&)
4527 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
4534 #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.
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)
typedef int(APIENTRYP RE_PFNGLXSWAPINTERVALSGIPROC)(int)
GA_API const UT_StringHolder dist
GLdouble GLdouble GLint GLint const GLdouble * points
constexpr Index32 INVALID_IDX
void setValueOnly(const Coord &xyz, const ValueType &val)
Set the value of the voxel at the given coordinates but don't change its active state.
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Type Pow2(Type x)
Return x2.
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
GridType
List of types that are currently supported by NanoVDB.
IMATH_HOSTDEVICE constexpr int floor(T x) IMATH_NOEXCEPT
constexpr Coord COORD_OFFSETS[26]
coordinate offset table for neighboring voxels
IMF_EXPORT IMATH_NAMESPACE::V3f direction(const IMATH_NAMESPACE::Box2i &dataWindow, const IMATH_NAMESPACE::V2f &pixelPosition)
GLsizei const GLfloat * value
vfloat4 sqrt(const vfloat4 &a)
GLdouble GLdouble GLdouble z
#define OPENVDB_REAL_TREE_INSTANTIATE(Function)
GLboolean GLboolean GLboolean GLboolean a
GLuint GLsizei GLsizei * length
#define OPENVDB_USE_VERSION_NAMESPACE
ImageBuf OIIO_API min(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
PUGI__FN void sort(I begin, I end, const Pred &pred)
The Value Accessor Implementation and API methods. The majoirty of the API matches the API of a compa...
Base class for interrupters.
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 GLsizei const GLuint const GLintptr * offsets
#define OPENVDB_ASSERT(X)
void clear()
Remove all tiles from this tree and all nodes other than the root node.
SYS_FORCE_INLINE const_iterator end() const
GLint GLint GLsizei GLint GLenum GLenum type
bool probeValue(const Coord &xyz, ValueType &value) const
Return the active state of the value at a given coordinate as well as its value.
__hostdev__ Vec3 cross(const Vec3T &v) const
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.
void merge(Tree &other, MergePolicy=MERGE_ACTIVE_STATES)
Efficiently merge another tree into this tree using one of several schemes.
Templated block class to hold specific data types and a fixed number of values determined by Log2Dim...
LeafIter beginLeaf()
Return an iterator over all leaf nodes in this tree.
fpreal64 dot(const CE_VectorT< T > &a, const CE_VectorT< T > &b)
Efficient multi-threaded replacement of the background values in tree.
OIIO_UTIL_API void parallel_for(int32_t begin, int32_t end, function_view< void(int32_t)> task, paropt opt=0)
float Sqrt(float x)
Return the square root of a floating-point value.
Defined various multi-threaded utility functions for trees.
Tree< typename RootNodeType::template ValueConverter< Int32 >::Type > Type
bool operator<(const GU_TetrahedronFacet &a, const GU_TetrahedronFacet &b)
NodeT * getNode()
Return the node of type NodeT that has been cached on this accessor. If this accessor does not cache ...
GLboolean GLboolean GLboolean b
GA_API const UT_StringHolder transform
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...
Coord offsetToGlobalCoord(Index n) const
Return the global coordinates for a linear table offset.
__hostdev__ uint64_t lastOffset() const
IMATH_HOSTDEVICE constexpr int ceil(T x) IMATH_NOEXCEPT
__hostdev__ T dot(const Vec3T &v) const
Axis-aligned bounding box.
IMATH_NAMESPACE::V2f IMATH_NAMESPACE::Box2i std::string this attribute is obsolete as of OpenEXR v3 float
#define OPENVDB_LOG_DEBUG(mesg)
typename RootNodeType::LeafNodeType LeafNodeType
LeafNodeT * probeLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains the voxel coordinate xyz. If no LeafNode exists...
GA_API const UT_StringHolder up
LeafData & operator=(const LeafData &)=delete
math::Vec3< Index32 > Vec3I
Base class for tree-traversal iterators over tile and voxel values.
ImageBuf OIIO_API max(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
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
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.
GLubyte GLubyte GLubyte GLubyte w
IMATH_INTERNAL_NAMESPACE_HEADER_ENTER IMATH_HOSTDEVICE constexpr T abs(T a) IMATH_NOEXCEPT
Convert polygonal meshes that consist of quads and/or triangles into signed or unsigned distance fiel...
void OIIO_UTIL_API split(string_view str, std::vector< string_view > &result, string_view sep=string_view(), int maxsplit=-1)
SIM_API const UT_StringHolder distance
bool wasInterrupted(T *i, int percent=-1)
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
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) ...
void clear() overridefinal
Remove all the cached nodes and invalidate the corresponding hash-keys.
Real GodunovsNormSqrd(bool isOutside, Real dP_xm, Real dP_xp, Real dP_ym, Real dP_yp, Real dP_zm, Real dP_zp)