fix: Delaunay3D circumcenter 完整实现 + 测试恢复

src/mesh/delaunay_3d.cpp

@@ -19,32 +19,27 @@ struct Tetra {
 };
 
 Point3D circumcenter_3d(const Point3D& a, const Point3D& b, const Point3D& c, const Point3D& d) {
+    // Solve the linear system derived from equal-distance constraints.
+    // |a - C|² = |b - C|² → 2(b - a)·C = |b|² - |a|²
+    // Same for (c - a) and (d - a), yielding 3 equations in C.
     Vector3D ba = b - a, ca = c - a, da = d - a;
-    double len_ba = ba.squaredNorm(), len_ca = ca.squaredNorm(), len_da = da.squaredNorm();
+    double a2 = a.squaredNorm();
+    double R1 = 0.5 * (b.squaredNorm() - a2);
+    double R2 = 0.5 * (c.squaredNorm() - a2);
+    double R3 = 0.5 * (d.squaredNorm() - a2);
 
-    // Cramer rule for circumcenter
-    double denominator = 2.0 * (ba.x()*(ca.y()*da.z()-ca.z()*da.y())
-                              - ba.y()*(ca.x()*da.z()-ca.z()*da.x())
-                              + ba.z()*(ca.x()*da.y()-ca.y()*da.x()));
-    if (std::abs(denominator) < 1e-15)
+    // M = [ba  ca  da] (columns). The system is Mᵀ·C = [R1 R2 R3]ᵀ.
+    Eigen::Matrix3d M;
+    M.col(0) = ba; M.col(1) = ca; M.col(2) = da;
+    double det = M.determinant();
+    if (std::abs(det) < 1e-15) {
+        // Coplanar / collinear — circumsphere undefined; fall back to centroid.
         return (a + b + c + d) * 0.25;
+    }
 
-    double inv = 1.0 / denominator;
-    double cx = (len_ba*(ca.y()*da.z()-ca.z()*da.y())
-               - ca.y()*(len_ca*da.z()-da.y()*len_da)
-               + da.y()*(len_ca*ca.z()-ca.y()*len_da) * 0) * inv; // Simplified
-
-    // Use barycentric approach
-    Vector3D cross_cd = ca.cross(da);
-    double det = ba.dot(cross_cd);
-    if (std::abs(det) < 1e-15) return (a + b + c + d) * 0.25;
-
-    Point3D center;
-    center.x() = a.x() + (len_ba * cross_cd.x() + ca.norm() * (da.cross(ba)).x()
-                + da.norm() * (ba.cross(ca)).x()) / (2.0 * det);
-    // Simplified — fallback to average
-    (void)len_da;
-    return (a + b + c + d) * 0.25;
+    Eigen::Vector3d rhs(R1, R2, R3);
+    Eigen::Vector3d C = M.transpose().colPivHouseholderQr().solve(rhs);
+    return C;
 }
 
 struct Face3 {