robust for collinear vertices of a face

algoim/organizer/organizer.hpp

@@ -344,7 +344,7 @@ void buildOcTreeV0(const Scene&             scene,
         int bbb = 1;
     }
     const std::vector<int>& polyIntersectIndices = node.polyIntersectIndices;
-    if (polyIntersectIndices.size() <= 8) {
+    if (polyIntersectIndices.size() <= 7) {
         leaves.emplace_back(node);
         return;
     }
@@ -545,7 +545,7 @@ BasicTaskRes basicTask(const Scene& scene, const OcTreeNode& node, int q = 10)
 
     uvector3i numSamples(50, 50, 50);
     // if (false) {
-    if (node.polyIntersectIndices.size() > 0) {
+    if (node.polyIntersectIndices.size() > 8) {
         // TODO: 蒙特卡洛
         std::cout << "Monte Carlo 1" << std::endl;
         ;

algoim/organizer/primitive.hpp

@@ -801,6 +801,28 @@ void makeHalfPlane(const HalfPlaneDesc& halfPlaneDesc, VisiblePrimitiveRep& visi
     buildNearBalancedBlobTree(visiblePrimitive.subBlobTree, 1);
 }
 
+uvector3 getFaceNorm(const std::vector<uvector3>& points, const std::vector<int>& indices, int faceBias, int faceStride)
+{
+    // 注意任意三点可能共线！
+
+    assert(faceStride >= 3);
+
+    // 遍历点集，寻找不共线的三个点
+    for (size_t i = 0; i < faceStride - 2; ++i) {
+        for (size_t j = i + 1; j < faceStride - 1; ++j) {
+            for (size_t k = j + 1; k < faceStride; ++k) {
+                uvector3 v1 = points[indices[faceBias + j]] - points[indices[faceBias + i]];
+                uvector3 v2 = points[indices[faceBias + k]] - points[indices[faceBias + i]];
+                uvector3 N  = cross(v1, v2);
+                if (!isnan(N) && norm(N) > std::numeric_limits<real>::epsilon()) {
+                    return N; // 返回归一化后的法向量
+                }
+            }
+        }
+    }
+    throw std::runtime_error("所有点都共线，无法计算法向量");
+}
+
 void makeMesh(const MeshDesc& mesh, VisiblePrimitiveRep& visiblePrimitive)
 {
     assert(visiblePrimitive.tensors.size() == mesh.indexInclusiveScan.size());
@@ -814,7 +836,7 @@ void makeMesh(const MeshDesc& mesh, VisiblePrimitiveRep& visiblePrimitive)
         auto&    indices                  = mesh.indices;
         uvector3 V01                      = vertices[indices[indexBeg + 1]] - mesh.vertices[indices[indexBeg]];
         uvector3 V02                      = vertices[indices[indexBeg + 2]] - mesh.vertices[indices[indexBeg]];
-        uvector3 N                        = cross(V01, V02);
+        uvector3 N                        = getFaceNorm(vertices, indices, indexBeg, indexSize);
         N                                /= norm(N);
         real d                            = -dot(N, vertices[indices[indexBeg]]);
         // 法线所指方向为>0区域

algoim/quadrature_multipoly.hpp

@@ -36,7 +36,7 @@
 #include "tanhsinh.hpp"
 #include "bernstein.hpp"
 
-#define STOP_WHEN_BLOCKED true
+#define STOP_WHEN_BLOCKED false
 
 #if STOP_WHEN_BLOCKED
 #include <chrono>

algoim/uvector.hpp

@@ -349,6 +349,14 @@ auto norm(const T& u)
     return sqrt(sqrnorm(u));
 }
 
+template <typename T, int N>
+bool isnan(const uvector<T, N>& u)
+{
+    for (int i = 0; i < N; ++i)
+        if (std::isnan(u(i))) return true;
+    return false;
+}
+
 // alter the extent of a uvector by truncation or zero backfill
 template <int M, typename T, int N>
 uvector<T, M> alter_extent(const uvector<T, N>& u)