subdivision

algoim/organizer/organizer.hpp

@@ -223,9 +223,9 @@ int binaryToDecimal(const uvector<int, N>& v, int zeroRep = 0)
 }
 
 template <int N>
-int replaceFirst(uvector<int, N>& v, int oldVal, int newVal)
+int replaceFirst(uvector<int, N>& 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<int, N>& v, int oldVal, int newVal)
     return -1;
 }
 
-void build(const Scene& scene, std::vector<Node>& nodes, int nowNodeIdx, const uvector3& nowNodeMin, const uvector3& nowNodeMax)
+template <int N>
+int findFirst(const uvector<int, N>& v, int val, int startIdx = 0)
+{
+    for (int i = startIdx; i < N; ++i) {
+        if (v(i) == val) return i;
+    }
+    return -1;
+}
+
+std::pair<uvector<int, 3>, uvector<int, 3>> getOneEightCellAABB(const uvector3&       min,
+                                                                const uvector3&       max,
+                                                                const uvector<int, 3> side,
+                                                                int                   negativeRep = -1)
+{
+    std::pair<uvector<int, 3>, uvector<int, 3>> 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<int, 2> 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<int, 3>&   mark,
+                                  std::vector<Node>& 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<Node>& nodes, int nowNodeIdx)
 {
     const std::vector<int>& 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<int, 3> 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) {
+            int signOnHalfPlane = bernstein::uniformSign(restrictToCenterPlane);
-                // primitive contain the conterface
+            if (signOnHalfPlane == -1) {
+                // primitive contain the centerface
                 mark(faceAxis) = 2;
-            } else {
+            } else if (signOnHalfPlane == 1) {
-                // primitive is on one side of the centerface
+                // 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());
+                tensor3 halfCellPoly(nullptr, poly.ext());
-                algoim_spark_alloc(real, negativeHalfCell);
+                algoim_spark_alloc(real, halfCellPoly);
-                bernstein::deCasteljau(poly, halfCellMin, halfCellMax, negativeHalfCell);
+                bernstein::deCasteljau(poly, halfCellMin, halfCellMax, halfCellPoly);
-                if (bernstein::uniformSign(negativeHalfCell) == -1) {
+                if (bernstein::uniformSign(halfCellPoly) != 1) {
-                    // 不在负半空间，此时一定在正半空间，可以assert一下正半空间的sign，不能为-1
+                    // 负空间有
-                    mark(faceAxis) = 1;
-                } else {
                     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<Node>& 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

algoim/organizer/primitive.hpp

@@ -595,11 +595,15 @@ struct Scene {
     // uvector3             min, max;
 };
 
+const int CHILD_NUM = 8;
+
 struct Node {
     std::vector<int> polyIntersectIndices;      // 同一poly会出现在不同node的polyIndices中
     std::vector<int> polyFullyContainedIndices; // 完全包含该node的poly, 同样会出现在不同node的polyIndices中
-    int              children[8];               // octree
+    std::array<int, CHILD_NUM> children;        // octree
-    uvector3         min, max;
+    uvector3                   min, max;
+
+    Node() { children.fill(-1); }
 };
 
 } // namespace algoim::Organizer