fix bug of blobtree traversal

algoim/organizer/blobtree.hpp

@@ -117,15 +117,15 @@ void propagate(BlobTree& tree, unsigned int nodeIdx, bool in)
     return;
 }
 
-// TODO: std::vector<bool> 可以改成bitset之类
+// TODO: std::vector<char> 浪费内存，应实现动态大小的bitset
 /**
  * relatedPrimitives 是当次traversal需要care的primitives的索引，例如OcTree子问题中关联的Primitive
  */
-int traverse(BlobTree& tree, const std::vector<int>& relatedPrimitives, const std::vector<bool>& primitiveInout)
+int traverse(BlobTree& tree, const std::vector<int>& relatedPrimitives, const std::vector<char>& primitiveInout)
 {
     assert(relatedPrimitives.size() == primitiveInout.size());
-    for (int primitiveIdx : relatedPrimitives) {
+    for (int i = 0; i < relatedPrimitives.size(); ++i) {
-        propagate(tree, tree.primitiveNodeIdx[primitiveIdx], primitiveInout[primitiveIdx]);
+        propagate(tree, tree.primitiveNodeIdx[relatedPrimitives[i]], static_cast<bool>(primitiveInout[i]));
         if (tree.structure.back().inOut != NODE_IN_OUT_UNKNOWN) { return tree.structure.back().inOut; }
     }
     return NODE_IN_OUT_UNKNOWN;

algoim/organizer/organizer.hpp

@@ -113,15 +113,20 @@ bool keepQuadraturePoint(const std::vector<tensor3>& tensors, const OcTreeNode&
     // 只需要考虑intersect polys，不用考虑fully contained polys
     const auto& polyIntersectIndices = ocTreeNode.polyIntersectIndices;
     assert(tensors.size() == polyIntersectIndices.size());
-    if (polyIntersectIndices.size() == 0) {
+    // if (polyIntersectIndices.size() == 0) {
-        assert(ocTreeNode.blobTree.structure.back().inOut != NODE_IN_OUT_UNKNOWN);
+    //     assert(ocTreeNode.blobTree.structure.back().inOut != NODE_IN_OUT_UNKNOWN);
-        return ocTreeNode.blobTree.structure.back().inOut;
+    //     return ocTreeNode.blobTree.structure.back().inOut == organizer::NODE_IN;
-    }
+    // }
-    std::vector<bool> primitiveInOuts(polyIntersectIndices.size());
+    std::vector<char> primitiveInOuts(polyIntersectIndices.size());
     for (int i = 0; i < ocTreeNode.polyIntersectIndices.size(); ++i) {
         // primitiveInOuts[i] = isInsidePowers(scene.polys[polyIntersectIndices[i]].rawPower, originPt);
         primitiveInOuts[i] = isInsideBernstein(tensors[i], point);
     }
+    if (primitiveInOuts.size() == 3 && primitiveInOuts[0] == NODE_IN && primitiveInOuts[1] == NODE_IN
+        && primitiveInOuts[2] == NODE_IN) {
+        int aaa = 1;
+        int bbb = 1;
+    }
     // 这里blobTree是拷贝传参
     auto blobTree = ocTreeNode.blobTree;
     int  res      = organizer::traverse(blobTree, ocTreeNode.polyIntersectIndices, primitiveInOuts);
@@ -615,15 +620,15 @@ void quadratureScene(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitiv
     for (int i = 0; i < minimalReps.size(); ++i) { rootNode.polyIntersectIndices.emplace_back(i); }
     int cnt = 1;
     buildOcTreeV0(scene, rootNode, leaves, 1, cnt);
-
+    std::cout << "octree built over" << std::endl;
-    basicTask(scene, leaves[14], q);
+    // basicTask(scene, leaves[14], q);
 
     int i = 0;
     for (const auto& leaf : leaves) {
         auto basicRes = basicTask(scene, leaf, q);
         if (std::isinf(basicRes.volume)) { std::cout << "inf volume when solving leaf: " << i << std::endl; }
         volume += basicRes.volume;
-        std::cout << "Solved leaves: " << ++i << "/" << leaves.size() << std::endl;
+        std::cout << "Solved leaves: " << ++i << "/" << leaves.size() << ", volume: " << volume << std::endl;
     }
 
     volume *= prod(xmax - xmin);

gjj/primitiveDebug.hpp

@@ -176,6 +176,23 @@ void caseScene()
     quadratureScene(primitiveDescriptions, uvector3(-0.8, -1.3, -1.6), uvector3(1.6, 1.6, 1.5), blobTree);
 }
 
+void caseScene1()
+{
+    const int                                   PRIMITIVE_CNT = 2;
+    std::vector<std::shared_ptr<PrimitiveDesc>> primitiveDescriptions(PRIMITIVE_CNT);
+    primitiveDescriptions[0] = std::make_shared<CuboidDesc>(CuboidDesc(0., 1.6));
+    primitiveDescriptions[1] = std::make_shared<CuboidDesc>(CuboidDesc(uvector3(0.6, 0.6, -0.6), 2.));
+    organizer::BlobTree blobTree;
+    blobTree.structure.resize(PRIMITIVE_CNT * 2 - 1);
+    blobTree.primitiveNodeIdx.resize(PRIMITIVE_CNT);
+    // blobTree.structure[0]     = {0, 0, NODE_IN_OUT_UNKNOWN, 0, 1, 2}; // cube1
+    blobTree.structure[0]     = {1, 0, NODE_IN_OUT_UNKNOWN, 0, 1, 2};               // cube1
+    blobTree.structure[1]     = {1, 0, NODE_IN_OUT_UNKNOWN, 0, 0, 0};               // cube2
+    blobTree.structure[2]     = {0, OP_INTERSECTION, NODE_IN_OUT_UNKNOWN, 0, 0, 0}; // INTERSECTION of cubes = opNode1
+    blobTree.primitiveNodeIdx = {0, 1};
+    quadratureScene(primitiveDescriptions, uvector3(-0.8, -1.3, -1.6), uvector3(1.6, 1.6, 1.5), blobTree);
+}
+
 void testDeCasteljau()
 {
     auto    phiDesc = std::make_shared<SphereDesc>(SphereDesc(0.8, uvector3(0), 1.));
@@ -343,7 +360,8 @@ void testPrimitive()
     // casePolyhedronSphere();
     // testSubDivideWithDeCasteljau();
     // testBlob();
-    caseScene();
+    // caseScene();
+    caseScene1();
     // testPlaneWithinSubcell();
     // caseCubeSphere();
     // caseTwoCube();