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