diff --git a/algoim/organizer/organizer.hpp b/algoim/organizer/organizer.hpp index d6ad4f6..e009084 100644 --- a/algoim/organizer/organizer.hpp +++ b/algoim/organizer/organizer.hpp @@ -223,9 +223,9 @@ int binaryToDecimal(const uvector& v, int zeroRep = 0) } template -int replaceFirst(uvector& v, int oldVal, int newVal) +int replaceFirst(uvector& v, int oldVal, int newVal, int startIdx = 0) { - for (int i = 0; i < N; ++i) { + for (int i = startIdx; i < N; ++i) { if (v(i) == oldVal) { v(i) = newVal; return i; @@ -234,42 +234,129 @@ int replaceFirst(uvector& v, int oldVal, int newVal) return -1; } -void build(const Scene& scene, std::vector& nodes, int nowNodeIdx, const uvector3& nowNodeMin, const uvector3& nowNodeMax) +template +int findFirst(const uvector& v, int val, int startIdx = 0) +{ + for (int i = startIdx; i < N; ++i) { + if (v(i) == val) return i; + } + return -1; +} + +std::pair, uvector> getOneEightCellAABB(const uvector3& min, + const uvector3& max, + const uvector side, + int negativeRep = -1) +{ + std::pair, uvector> res = {min, max}; + uvector3 mid = (min + max) / 2; + for (int i = 0; i < 3; ++i) { + if (side(i) == negativeRep) { + res.second(i) = mid(i); + } else { + res.first(i) = mid(i); + } + } + return res; +} + +std::array sides = {-1, 1}; + +// 对于mark含2个及以上0的情况，尝试对每个0填1或-1的所有组合 +// TODO: 参数太多了，考虑换用std::function + lambda +void visitSubcellOnBothSidesOfDir(const uvector3& nowNodeMin, + const uvector3& nowNodeMax, + const uvector3& nodeMid, + const int lastIdx, + const int polyIntersectIndex, + const tensor3& poly, + + int startIdxToCheck, + uvector& mark, + std::vector& nodes) +{ + int zeroIdx = findFirst(mark, 0, startIdxToCheck); + if (zeroIdx == -1) { + tensor3 halfCellPoly(nullptr, poly.ext()); + algoim_spark_alloc(real, halfCellPoly); + auto subNodeMin = nowNodeMin; + auto subNodeMax = nowNodeMax; + for (int k = 0; k < 3; ++k) { + if (mark(k) == -1) { + subNodeMax(k) = nodeMid(k); + } else if (mark(k) == 1) { + subNodeMin(k) = nodeMid(k); + } else { + std::cerr << "error: mark = " << mark; + } + } + bernstein::deCasteljau(poly, subNodeMin, subNodeMax, halfCellPoly); // 求1/8空间下的表达 + if (bernstein::uniformSign(halfCellPoly) != 1) { + nodes[lastIdx + binaryToDecimal(mark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndex); + } + } else { + for (auto side : sides) { + mark(zeroIdx) = side; + visitSubcellOnBothSidesOfDir(nowNodeMin, + nowNodeMax, + nodeMid, + lastIdx, + polyIntersectIndex, + poly, + zeroIdx + 1, + mark, + nodes); + } + } +} + +void build(const Scene& scene, std::vector& nodes, int nowNodeIdx) { const std::vector& polyIntersectIndices = nodes[nowNodeIdx].polyIntersectIndices; if (polyIntersectIndices.size() == 1) { // TODO 相交很少时 } - int lastIdx = nodes.size(); + const uvector3& nowNodeMin = nodes[nowNodeIdx].min; + const uvector3& nowNodeMax = nodes[nowNodeIdx].max; + int lastIdx = nodes.size(); nodes.resize(lastIdx + 8); uvector3 nodeMid = (nowNodeMin + nowNodeMax) / 2; for (int i = 0; i < polyIntersectIndices.size(); ++i) { - auto poly = scene.polys[polyIntersectIndices[i]]; + const auto& poly = scene.polys[polyIntersectIndices[i]]; uvector mark(0, 0, 0); for (int faceAxis = 0; faceAxis < 3; ++faceAxis) { real centerPlane = nodeMid(faceAxis); tensor2 restrictToCenterPlane(nullptr, remove_component(poly.ext(), faceAxis)); algoim_spark_alloc(real, restrictToCenterPlane); restrictToInnerFace(poly, faceAxis, centerPlane, restrictToCenterPlane); - if (bernstein::uniformSign(restrictToCenterPlane) == -1) { - // primitive contain the conterface + int signOnHalfPlane = bernstein::uniformSign(restrictToCenterPlane); + if (signOnHalfPlane == -1) { + // primitive contain the centerface mark(faceAxis) = 2; - } else { - // primitive is on one side of the centerface + } else if (signOnHalfPlane == 1) { + // primitive intersects either side or both sides of the centerface // deCasteljau to transformation + // TODO: not 0, 1, or 0.5 uvector3 halfCellMin = 0, halfCellMax = 1; halfCellMax(faceAxis) = 0.5; - tensor3 negativeHalfCell(nullptr, poly.ext()); - algoim_spark_alloc(real, negativeHalfCell); - bernstein::deCasteljau(poly, halfCellMin, halfCellMax, negativeHalfCell); - if (bernstein::uniformSign(negativeHalfCell) == -1) { - // 不在负半空间，此时一定在正半空间，可以assert一下正半空间的sign，不能为-1 - mark(faceAxis) = 1; - } else { + tensor3 halfCellPoly(nullptr, poly.ext()); + algoim_spark_alloc(real, halfCellPoly); + bernstein::deCasteljau(poly, halfCellMin, halfCellMax, halfCellPoly); + if (bernstein::uniformSign(halfCellPoly) != 1) { + // 负空间有 mark(faceAxis) = -1; } + // TODO: not 0, 1, or 0.5 + halfCellMax(faceAxis) = 1; + halfCellMin(faceAxis) = 0.5; + bernstein::deCasteljau(poly, halfCellMin, halfCellMax, halfCellPoly); + if (bernstein::uniformSign(halfCellPoly) != 1) { + // 正空间有 + mark(faceAxis) += 1; // 当正负空间都有，记0 + } } } + if (any(mark == uvector3(2, 2, 2))) { int subIdx = 0; for (MultiLoop<3> j(0, 2); ~j; ++j, ++subIdx) { @@ -298,8 +385,28 @@ void build(const Scene& scene, std::vector& nodes, int nowNodeIdx, const u nodes[lastIdx + binaryToDecimal(sideMark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndices[i]); } else { // zeroCount == 2 or 3 + visitSubcellOnBothSidesOfDir(nowNodeMin, + nowNodeMax, + nodeMid, + lastIdx, + polyIntersectIndices[i], + poly, + 0, + mark, + nodes); } } + + // launch subdivision in subcells + int subIdx = 0; + for (MultiLoop<3> j(0, 2); ~j; ++j, ++subIdx) { + int idx = lastIdx + subIdx; + auto nodeAABB = getOneEightCellAABB(nodes[lastIdx + subIdx].min, nodes[lastIdx + subIdx].max, j(), 0); + nodes[idx].min = nodeAABB.first; + nodes[idx].max = nodeAABB.second; + nodes[nowNodeIdx].children[subIdx] = idx; + build(scene, nodes, idx); + } // TODO: 考虑fully contain 问题 } }; // namespace algoim::Organizer diff --git a/algoim/organizer/primitive.hpp b/algoim/organizer/primitive.hpp index 66aca3c..e218387 100644 --- a/algoim/organizer/primitive.hpp +++ b/algoim/organizer/primitive.hpp @@ -595,11 +595,15 @@ struct Scene { // uvector3 min, max; }; +const int CHILD_NUM = 8; + struct Node { std::vector polyIntersectIndices; // 同一poly会出现在不同node的polyIndices中 std::vector polyFullyContainedIndices; // 完全包含该node的poly, 同样会出现在不同node的polyIndices中 - int children[8]; // octree - uvector3 min, max; + std::array children; // octree + uvector3 min, max; + + Node() { children.fill(-1); } }; } // namespace algoim::Organizer \ No newline at end of file