From f44e3efe477b6d7af9064e013ad7eb64496537ba Mon Sep 17 00:00:00 2001 From: gjj Date: Sun, 1 Sep 2024 23:32:46 +0800 Subject: [PATCH] subdivision --- algoim/organizer/organizer.hpp | 123 ++++++++++++--------------------- algoim/organizer/primitive.hpp | 18 +++-- gjj/primitiveDebug.hpp | 32 ++++++++- 3 files changed, 89 insertions(+), 84 deletions(-) diff --git a/algoim/organizer/organizer.hpp b/algoim/organizer/organizer.hpp index e009084..e63fcc7 100644 --- a/algoim/organizer/organizer.hpp +++ b/algoim/organizer/organizer.hpp @@ -264,65 +264,55 @@ 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, +void visitSubcellOnBothSidesOfDir(const uvector3& nodeMid, const int polyIntersectIndex, const tensor3& poly, - int startIdxToCheck, - uvector& mark, - std::vector& nodes) + std::array& subNodes, + int startIdxToCheck, + uvector& mark) { 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); - } + int subIdx = binaryToDecimal(mark); + auto& subNode = subNodes[subIdx]; + bernstein::deCasteljau(poly, subNode.min, subNode.max, halfCellPoly); // 求1/8空间下的表达 + if (bernstein::uniformSign(halfCellPoly) != 1) { subNode.polyIntersectIndices.emplace_back(polyIntersectIndex); } } else { for (auto side : sides) { mark(zeroIdx) = side; - visitSubcellOnBothSidesOfDir(nowNodeMin, - nowNodeMax, - nodeMid, - lastIdx, - polyIntersectIndex, - poly, - zeroIdx + 1, - mark, - nodes); + visitSubcellOnBothSidesOfDir(nodeMid, polyIntersectIndex, poly, subNodes, zeroIdx + 1, mark); } + mark(zeroIdx) = 0; } } -void build(const Scene& scene, std::vector& nodes, int nowNodeIdx) +// 叶节点单独存在leaves中，是因为最后只需要叶节点信息，没有自顶向下/自底向上遍历 +// 所以树结构是不需要的记录的 +// void build(const Scene& scene, std::vector& intermediateNodes, std::vector leaves, int nowNodeIdx) +void build(const Scene& scene, const Node& node, std::vector leaves) { - const std::vector& polyIntersectIndices = nodes[nowNodeIdx].polyIntersectIndices; - if (polyIntersectIndices.size() == 1) { - // TODO 相交很少时 + const std::vector& polyIntersectIndices = node.polyIntersectIndices; + if (polyIntersectIndices.size() <= 4) { + leaves.emplace_back(node); + return; + } + const uvector3& nowNodeMin = node.min; + const uvector3& nowNodeMax = node.max; + std::array subNodes; + // intermediateNodes.resize(lastIdx + 8); + int subIdx = 0; + for (MultiLoop<3> j(0, 2); ~j; ++j, ++subIdx) { + auto nodeAABB = getOneEightCellAABB(node.min, node.max, j(), 0); + subNodes[subIdx].min = nodeAABB.first; + subNodes[subIdx].max = nodeAABB.second; } - 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) { - const auto& poly = scene.polys[polyIntersectIndices[i]]; + const int polyIntersectIndex = polyIntersectIndices[i]; + const auto& poly = scene.polys[polyIntersectIndex]; uvector mark(0, 0, 0); for (int faceAxis = 0; faceAxis < 3; ++faceAxis) { real centerPlane = nodeMid(faceAxis); @@ -336,9 +326,8 @@ void build(const Scene& scene, std::vector& nodes, int nowNodeIdx) } 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; + uvector3 halfCellMin = nowNodeMin, halfCellMax = nowNodeMax; + halfCellMax(faceAxis) = nodeMid(faceAxis); tensor3 halfCellPoly(nullptr, poly.ext()); algoim_spark_alloc(real, halfCellPoly); bernstein::deCasteljau(poly, halfCellMin, halfCellMax, halfCellPoly); @@ -346,9 +335,8 @@ void build(const Scene& scene, std::vector& nodes, int nowNodeIdx) // 负空间有 mark(faceAxis) = -1; } - // TODO: not 0, 1, or 0.5 - halfCellMax(faceAxis) = 1; - halfCellMin(faceAxis) = 0.5; + halfCellMax(faceAxis) = nowNodeMax(faceAxis); + halfCellMin(faceAxis) = nodeMid(faceAxis); bernstein::deCasteljau(poly, halfCellMin, halfCellMax, halfCellPoly); if (bernstein::uniformSign(halfCellPoly) != 1) { // 正空间有 @@ -358,55 +346,32 @@ void build(const Scene& scene, std::vector& nodes, int nowNodeIdx) } if (any(mark == uvector3(2, 2, 2))) { - int subIdx = 0; - for (MultiLoop<3> j(0, 2); ~j; ++j, ++subIdx) { - nodes[lastIdx + subIdx].polyIntersectIndices.emplace_back(poly); - nodes[lastIdx + subIdx].min = nowNodeMin; - nodes[lastIdx + subIdx].max = nowNodeMax; - for (int k = 0; k < 3; ++k) { - if (j(k) == 0) - nodes[lastIdx + subIdx].max(k) = nodeMid(k); - else - nodes[lastIdx + subIdx].min(k) = nodeMid(k); - } + for (int subIdx = 0; subIdx < CHILD_NUM; ++subIdx) { + // intermediateNodes[lastIdx + subIdx].polyIntersectIndices.emplace_back(poly); + // intermediateNodes[lastIdx + subIdx].min = nowNodeMin; + // intermediateNodes[lastIdx + subIdx].max = nowNodeMax; + subNodes[subIdx].polyIntersectIndices.emplace_back(polyIntersectIndex); } continue; } int zeroCount = numOfZero<3>(mark); if (zeroCount == 0) { // mark has -1 or 1 only - nodes[lastIdx + binaryToDecimal(mark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndices[i]); + subNodes[binaryToDecimal(mark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndex); } else if (zeroCount == 1) { // poly related to 2 subcells uvector sideMark = mark; int zeroIdx = replaceFirst(sideMark, 0, 1); - nodes[lastIdx + binaryToDecimal(sideMark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndices[i]); + subNodes[binaryToDecimal(sideMark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndex); sideMark(zeroIdx) = -1; - nodes[lastIdx + binaryToDecimal(sideMark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndices[i]); + subNodes[binaryToDecimal(sideMark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndex); } else { // zeroCount == 2 or 3 - visitSubcellOnBothSidesOfDir(nowNodeMin, - nowNodeMax, - nodeMid, - lastIdx, - polyIntersectIndices[i], - poly, - 0, - mark, - nodes); + visitSubcellOnBothSidesOfDir(nodeMid, polyIntersectIndex, poly, subNodes, 0, mark); } } - // 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); - } + for (subIdx = 0; subIdx < CHILD_NUM; ++subIdx) { build(scene, subNodes[subIdx], leaves); } // TODO: 考虑fully contain 问题 } }; // namespace algoim::Organizer diff --git a/algoim/organizer/primitive.hpp b/algoim/organizer/primitive.hpp index e218387..494eb68 100644 --- a/algoim/organizer/primitive.hpp +++ b/algoim/organizer/primitive.hpp @@ -600,10 +600,20 @@ const int CHILD_NUM = 8; struct Node { std::vector polyIntersectIndices; // 同一poly会出现在不同node的polyIndices中 std::vector polyFullyContainedIndices; // 完全包含该node的poly, 同样会出现在不同node的polyIndices中 - std::array children; // octree - uvector3 min, max; - - Node() { children.fill(-1); } + // std::array children; // octree + uvector3 min, max; + + // Node() { children.fill(-1); } + Node() = default; + + // 深拷贝 + Node(const Node& node) + : polyIntersectIndices(node.polyIntersectIndices), + polyFullyContainedIndices(node.polyFullyContainedIndices), + min(node.min), + max(node.max) + { + } }; } // namespace algoim::Organizer \ No newline at end of file diff --git a/gjj/primitiveDebug.hpp b/gjj/primitiveDebug.hpp index cc30cdb..8962329 100644 --- a/gjj/primitiveDebug.hpp +++ b/gjj/primitiveDebug.hpp @@ -158,8 +158,38 @@ void casePolyhedronSphere() basicTask(primitiveDescriptions); } +void testDeCasteljau() +{ + auto phiDesc = std::make_shared(SphereDesc(0.8, uvector3(0), 1.)); + tensor3 tensor(nullptr, 3), tensor01(nullptr, 3), tensorBernstein(nullptr, 3), transformedTensorBernstein(nullptr, 3); + algoim_spark_alloc(real, tensor, tensor01, tensorBernstein, transformedTensorBernstein); + makeSphere(*phiDesc, tensor); + uvector3 xmax = 1, xmin = -1, range = xmax - xmin; + Organizer::detail::powerTransformation(range, xmin, tensor, tensor01); + Organizer::detail::power2BernsteinTensor(tensor01, tensorBernstein); + + // bernstein::deCasteljau(tensorBernstein, uvector3(0), uvector3(0.2690598923241497 + 0.0000001), + // transformedTensorBernstein); + bernstein::deCasteljau(tensorBernstein, + uvector3(0.2690598923241497 + 0.0000001), + uvector3(1 - 0.2690598923241497 - 0.0000001), + transformedTensorBernstein); + uvector3 testX = uvector3(0.5); + + int sign = bernstein::uniformSign(transformedTensorBernstein); + real evalX = evalPower(tensor, testX); + + std::cout << "evalX original rep = " << evalX << std::endl; + evalX = bernstein::evalBernsteinPoly(tensorBernstein, testX); + std::cout << "evalX bernstein within 0-1 = " << evalX << std::endl; + evalX = bernstein::evalBernsteinPoly(transformedTensorBernstein, testX); + std::cout << "evalX bernstein after Decasteljau = " << evalX << std::endl; + std::cout << "sign = " << sign << std::endl; +} + void testPrimitive() { - casePolyhedron1(); + // casePolyhedron1(); // casePolyhedronSphere(); + testDeCasteljau(); } \ No newline at end of file