使用CGAL计算SurfaceMesh主曲率及解决代码报错问题

一、先解决当前代码的异常问题 #

你遇到的nb_input_pts >= nb_d_jet_coeff异常，核心是输入拟合的点数量不足，加上参数传递错误，具体修正步骤：

1. 修正monge_fit的维度参数
   你调用monge_fit(adj.begin(), adj.end(), 2, 1)时，第四个参数是空间维度，3D网格场景下应该传3，不是1。这个参数决定了拟合的空间维度，传1会导致计算需要的系数数量错误，同时也不符合3D点的拟合场景。

2. 补充足够的邻域点

   • 首先要把当前顶点本身的点加入邻域列表：你的代码只收集了相邻顶点，没有包含顶点v自己的坐标，这会减少可用的拟合点数量。
   • 仅取1-ring相邻顶点可能不够（比如边界顶点的相邻顶点数量少于拟合所需），建议扩展到2-ring邻域，或者使用网格的面来收集邻域内的点（比如遍历顶点周围的面，把面的所有顶点都加入邻域，去重后使用）。

   修正后的邻域收集代码示例：

   adj.push_back(result.point(v)); // 先加入当前顶点
   auto h = result.halfedge(v);
   for (HalfedgeDescriptor he : result.halfedges_around_target(h))
   {
       Point p = result.point(result.source(he));
       adj.push_back(p);
       // 扩展到2-ring：遍历相邻顶点的邻域
       HalfedgeDescriptor he2 = result.halfedge(result.source(he));
       for (HalfedgeDescriptor he_around : result.halfedges_around_target(he2))
       {
           Point p2 = result.point(result.source(he_around));
           adj.push_back(p2);
       }
   }
   // 去重，避免重复点影响拟合
   std::sort(adj.begin(), adj.end());
   adj.erase(std::unique(adj.begin(), adj.end()), adj.end());
   

3. 增加点数量检查
   在调用monge_fit前，先判断邻域点数量是否满足要求，避免异常：

   // 对于jet order=2，3D下需要至少7个点（二次多项式有7个系数）
   if (adj.size() < 7) {
       // 跳过或用默认值处理，比如边界顶点
       continue;
   }
   

二、关于comply_wrt_given_normal的作用 #

这个函数是用来统一曲率的符号方向的：Jet拟合得到的Monge曲面的法线方向可能和你传入的顶点法线（vn.first[v]）相反，导致计算出的曲率符号颠倒（比如凸面算出负曲率）。调用这个函数后，Monge形式会调整自身的朝向，使得曲面法线和给定的顶点法线方向一致，保证曲率符号符合网格的几何直观。

三、更简便的批量曲率计算方法 #

CGAL提供了专门的顶点曲率估计函数，无需手动实现邻域收集和拟合，直接批量计算所有顶点的曲率：
使用CGAL::estimate_vertex_curvatures函数（需要包含头文件#include <CGAL/Surface_mesh_curvature_estimate.h>），示例代码：

#include <CGAL/Surface_mesh_curvature_estimate.h>

// 假设result是你的SurfaceMesh，Kernel是你使用的CGAL内核
std::map<Vertex_descriptor, CGAL::Curvature> curvatures;
CGAL::estimate_vertex_curvatures(result, std::inserter(curvatures, curvatures.begin()));

// 遍历获取每个顶点的曲率
for (auto v : result.vertices()) {
    double k1 = curvatures[v].max_principal_curvature;
    double k2 = curvatures[v].min_principal_curvature;
    double gaussian = curvatures[v].gaussian_curvature;
    double mean = curvatures[v].mean_curvature;
}

这个函数内部已经封装了邻域选择、Jet拟合等逻辑，能直接输出主曲率、高斯曲率、平均曲率，比手动实现更高效可靠。

内容的提问来源于stack exchange，提问作者senseiwa