From 028b0c3bc896ee7f8b50acfa25ad1f8885f999da Mon Sep 17 00:00:00 2001 From: gjj Date: Sun, 22 Sep 2024 17:35:27 +0800 Subject: [PATCH] fix bug of octree construction and blobtree propagation --- algoim/organizer/blobtree.hpp | 6 +++- algoim/organizer/organizer.hpp | 55 +++++++++++++++++++++++---------- algoim/organizer/primitive.hpp | 9 ++++-- algoim/quadrature_multipoly.hpp | 3 ++ gjj/primitiveDebug.hpp | 4 +-- 5 files changed, 56 insertions(+), 21 deletions(-) diff --git a/algoim/organizer/blobtree.hpp b/algoim/organizer/blobtree.hpp index d9db81a..f0c09c8 100644 --- a/algoim/organizer/blobtree.hpp +++ b/algoim/organizer/blobtree.hpp @@ -65,6 +65,10 @@ void propagate(BlobTree& tree, unsigned int nodeIdx, bool in) { const std::size_t rootIdx = tree.structure.size() - 1; auto& node = tree.structure[nodeIdx]; + if (nodeIdx == 4) { + int aaa = 1; + int bbb = 1; + } if (node.inOut != NODE_IN_OUT_UNKNOWN) { int aaa = 1; int bb = 1; @@ -73,7 +77,7 @@ void propagate(BlobTree& tree, unsigned int nodeIdx, bool in) node.inOut = in ? NODE_IN : NODE_OUT; for (unsigned int nowIdx = nodeIdx; nowIdx != rootIdx;) { const auto& nowNode = tree.structure[nowIdx]; - unsigned int nextIdx = isLeft(nowNode) ? nodeIdx + nowNode.ancestor : nodeIdx + 1; + unsigned int nextIdx = isLeft(nowNode) ? nowIdx + nowNode.ancestor : nowIdx + 1; auto& nextNode = tree.structure[nextIdx]; // father if (nextNode.inOut != 0) { return; } if (nowNode.inOut == NODE_OUT) { diff --git a/algoim/organizer/organizer.hpp b/algoim/organizer/organizer.hpp index 706ec49..f70912c 100644 --- a/algoim/organizer/organizer.hpp +++ b/algoim/organizer/organizer.hpp @@ -108,11 +108,7 @@ AABB getOneEightCellAABB(const AABB& fatherAABB, const uvector side, int } } // namespace detail -// 这里blobTree是拷贝传参 -bool keepQuadraturePoint(const Scene& scene, - organizer::BlobTree blobTree, - const OcTreeNode& ocTreeNode, - const uvector3& originPt) +bool keepQuadraturePoint(const Scene& scene, const OcTreeNode& ocTreeNode, const uvector3& originPt) { // 只需要考虑intersect polys，不用考虑fully contained polys const auto& polyIntersectIndices = ocTreeNode.polyIntersectIndices; @@ -121,7 +117,9 @@ bool keepQuadraturePoint(const Scene& scene, // primitiveInOuts[i] = isInsidePowers(scene.polys[polyIntersectIndices[i]].rawPower, originPt); primitiveInOuts[i] = isInsideBernstein(scene.tensors[polyIntersectIndices[i]], originPt); } - int res = organizer::traverse(blobTree, ocTreeNode.polyIntersectIndices, primitiveInOuts); + // 这里blobTree是拷贝传参 + auto blobTree = ocTreeNode.blobTree; + int res = organizer::traverse(blobTree, ocTreeNode.polyIntersectIndices, primitiveInOuts); assert(res == organizer::NODE_IN || res == organizer::NODE_OUT); return res == organizer::NODE_IN; } @@ -285,13 +283,13 @@ void buildOcTreeV1(const Scene& scene, const OcTreeNode& node, std::vector& leaves, int depth) +void buildOcTreeV0(const Scene& scene, const OcTreeNode& node, std::vector& leaves, int depth, int& cnt) { if (deepest < depth) { deepest = depth; - std::cout << "depth: " << depth << std::endl; + std::cout << "depth: " << deepest << std::endl; } if (depth == 14300) { int aaa = 1; @@ -302,6 +300,11 @@ void buildOcTreeV0(const Scene& scene, const OcTreeNode& node, std::vector subNodes; std::vector subNodes(CHILD_NUM); // intermediateNodes.resize(lastIdx + 8); @@ -311,6 +314,15 @@ void buildOcTreeV0(const Scene& scene, const OcTreeNode& node, std::vector(); } + if (node.blobTree.structure[4].inOut != NODE_IN_OUT_UNKNOWN + && std::find(polyIntersectIndices.begin(), polyIntersectIndices.end(), 3) != polyIntersectIndices.end()) { + int aaa = 1; + int bbb = 1; + } + if (cnt == 1) { + int aaa = 1; + int bbb = 1; + } uvector3 nodeMid = node.aabb.center(); for (int i = 0; i < polyIntersectIndices.size(); ++i) { const int polyIntersectIndex = polyIntersectIndices[i]; @@ -322,6 +334,10 @@ void buildOcTreeV0(const Scene& scene, const OcTreeNode& node, std::vector>& primitives, in // for (auto& p : originTensorStacks) delete p; }; -BasicTaskRes basicTask(const Scene& scene, const organizer::BlobTree& blobTree, const OcTreeNode& node, int q = 10) +BasicTaskRes basicTask(const Scene& scene, const OcTreeNode& node, int q = 10) { auto integrand = [](const uvector& x) { return 1.0; }; real volume = 0., surf = 0.; @@ -442,10 +464,10 @@ BasicTaskRes basicTask(const Scene& scene, const organizer::BlobTree& blobTree, ImplicitPolyQuadrature<3> ipquad(scene.tensors, node.polyIntersectIndices); ipquad.integrate(AutoMixed, q, [&](const uvector& x, real w) { auto realX = x * range + node.aabb.min; - if (keepQuadraturePoint(scene, blobTree, node, realX)) { volume += w * integrand(realX); } + if (keepQuadraturePoint(scene, node, realX)) { volume += w * integrand(realX); } }); ipquad.integrate_surf(AutoMixed, q, [&](const uvector& x, real w, const uvector& wn) { - surf += w * integrand(x * range + node.aabb.min); + volume += w * integrand(x * range + node.aabb.min); }); return {volume, surf}; } @@ -479,7 +501,7 @@ void quadratureScene(const std::vector>& primitiv } else if (auto pt = std::dynamic_pointer_cast(primitives[i])) { const int faceCount = pt->indexInclusiveScan.size(); assert(faceCount > 1); - visiblePrimitiveReps[i].tensors = std::vector(faceCount, tensor3(nullptr, 2)); + visiblePrimitiveReps[i].tensors = std::vector(faceCount, tensor3(nullptr, 3)); auto& planeTensors = visiblePrimitiveReps[i].tensors; algoimSparkAllocHeapVector(tensorStacks, planeTensors); @@ -566,10 +588,11 @@ void quadratureScene(const std::vector>& primitiv Scene scene{realPrimitives, AABB(xmin, xmax)}; OcTreeNode rootNode(0, 1, blobTree); for (int i = 0; i < realPrimitives.size(); ++i) { rootNode.polyIntersectIndices.emplace_back(i); } - buildOcTreeV0(scene, rootNode, leaves, 1); + int cnt = 1; + buildOcTreeV0(scene, rootNode, leaves, 1, cnt); for (const auto& leaf : leaves) { - auto basicRes = basicTask(scene, blobTree, leaf, q); + auto basicRes = basicTask(scene, leaf, q); volume += basicRes.volume * prod(leaf.aabb.max - leaf.aabb.min); } diff --git a/algoim/organizer/primitive.hpp b/algoim/organizer/primitive.hpp index 60743c8..8fb0f5d 100644 --- a/algoim/organizer/primitive.hpp +++ b/algoim/organizer/primitive.hpp @@ -699,12 +699,17 @@ struct VisiblePrimitiveRep { organizer::BlobTree subBlobTree; }; +struct MinimalPrimitiveRep { + tensor3 tensor; + AABB aabb; +}; + struct Scene { // std::vector polys; // std::vector visiblePrimitives; // organizer::BlobTree blobTree; - std::vector tensors; - AABB boundary; + std::vector tensors; + AABB boundary; }; void upwardGeneratingNodes(std::vector& levelLeftBook, diff --git a/algoim/quadrature_multipoly.hpp b/algoim/quadrature_multipoly.hpp index d199f4d..0979334 100644 --- a/algoim/quadrature_multipoly.hpp +++ b/algoim/quadrature_multipoly.hpp @@ -281,6 +281,9 @@ bool resultant_core(const xarray& p, const xarray* q, int k, x xarray mat(nullptr, uvector{M, M}); /**gjj***/ + // 当Q == 0，应该是对p求P类点，且p完全垂直降维方向的情形，这样的R类点可以直接丢弃 + // 因为产生的R类点和Q类点是重合的（平行降维方向的多项式一定和上下边界交于同一位置，其降维后也是R类点降维后的结果） + // 当P和Q都<=1，要么是Q == 0，要么是两个垂直降维方向的多项式求R类点，此时交线退化为点，也丢弃 if (Q == 0 || (P <= 1 && Q <= 1)) { return false; } if (!(P >= 2 && Q >= 1 && M >= 1)) { int aaa = 1; diff --git a/gjj/primitiveDebug.hpp b/gjj/primitiveDebug.hpp index 2339e27..2a3723b 100644 --- a/gjj/primitiveDebug.hpp +++ b/gjj/primitiveDebug.hpp @@ -140,8 +140,8 @@ void caseScene() const int PRIMITIVE_CNT = 6; std::vector> primitiveDescriptions(PRIMITIVE_CNT); primitiveDescriptions[0] = std::make_shared(CuboidDesc(0., 1.6)); - // primitiveDescriptions[1] = std::make_shared(CuboidDesc(uvector3(0.6, 0.6, -0.6), 2.)); - primitiveDescriptions[1] = std::make_shared(SphereDesc(1., uvector3(0.6, 0.6, -0.6), 1.)); + primitiveDescriptions[1] = std::make_shared(CuboidDesc(uvector3(0.6, 0.6, -0.6), 2.)); + // primitiveDescriptions[1] = std::make_shared(SphereDesc(1., uvector3(0.6, 0.6, -0.6), 1.)); primitiveDescriptions[2] = std::make_shared(SphereDesc(0.7, uvector3(0.8, 0.8, 0.8), 1.)); primitiveDescriptions[3] = std::make_shared(SphereDesc(0.5, uvector3(-0.3, -0.8, 0.8), 1.)); std::vector pyramidBottomVertices = {