基于CGAL实现变度量场下Delaunay三角剖分及高阶嵌入问询

一、变度量场下的Delaunay三角剖分（CGAL实现） #

CGAL没有直接提供Riemannian度量下的Delaunay剖分接口，但可以通过自定义几何traits类适配其三角剖分框架实现，核心是重定义Delaunay剖分依赖的几何判断逻辑：

• 核心思路：CGAL的Triangulation_2/Triangulation_3是基于traits类的通用框架，默认traits采用欧氏kernel，你可以实现自定义traits，将所有几何操作替换为Riemannian度量下的计算。
• 关键函数重写：
  • 必须实现side_of_oriented_circle：这个函数决定Delaunay剖分的核心条件——判断点是否在三角形的“外接圆”（Riemannian度量下的测地圆）内部/外部/边界。需要基于每个顶点关联的Riemannian度量矩阵，计算对应的测地圆方程，再判断点的位置。
  • 辅助函数：比如compute_squared_distance需要返回Riemannian度量下的平方距离，而非欧氏距离。
• 度量场的关联：将每个点对应的Riemannian度量（比如2x2/3x3对称正定矩阵）存储在可被traits访问的容器中（比如std::vector，按顶点索引对应），traits在计算时根据点的索引获取对应的度量参数。

示例伪代码逻辑：

// 自定义Riemannian度量的traits类
struct Riemannian_Traits : public CGAL::Triangulation_traits_2<CGAL::Exact_predicates_inexact_constructions_kernel> {
  // 重写side_of_oriented_circle函数
  template<typename Point, typename Circle>
  CGAL::Oriented_side side_of_oriented_circle(const Point& p, const Circle& c) const {
    // 获取p和c三个顶点对应的Riemannian度量矩阵
    const auto& metric_p = get_metric(p);
    const auto& metric_a = get_metric(c.point(0));
    const auto& metric_b = get_metric(c.point(1));
    const auto& metric_c = get_metric(c.point(2));
    // 计算Riemannian度量下的测地圆判断逻辑，返回对应的Oriented_side
    return compute_riemannian_side(metric_p, metric_a, metric_b, metric_c, p, c);
  }
};

// 初始化三角剖分
Riemannian_Traits traits;
CGAL::Triangulation_2<Riemannian_Traits> tri(traits);
// 插入带度量的点（需要自己维护点到度量的映射）
for (const auto& [point, metric] : point_metric_pairs) {
  auto vh = tri.insert(point);
  // 将metric存储到对应顶点的属性中
  vh->set_metric(metric);
}

二、给定网格与顶点度量场的高阶嵌入 #

针对你提到的BAMG网格丢失度量场的问题，以及高阶嵌入需求，可按以下方向处理：

1. 保留顶点度量场 #

不管用BAMG还是CGAL生成网格，都可以手动维护顶点与度量场的关联：

• 若从BAMG导出网格，可将顶点ID和对应的度量矩阵存储为附加数据（比如在OBJ文件中添加自定义属性行，或用JSON/CSV单独存储映射关系）；
• 导入CGAL后，使用CGAL::Surface_mesh的顶点属性机制，为每个顶点添加度量属性：

  CGAL::Surface_mesh<Point> mesh;
  // 导入网格后，添加顶点属性
  auto metric_prop = mesh.add_property_map<CGAL::Vertex_index, Eigen::Matrix2d>("v:metric").first;
  // 遍历顶点，关联度量场
  for (auto v : mesh.vertices()) {
    metric_prop[v] = get_metric_by_vertex_id(v.id());
  }
  

2. 高阶嵌入实现 #

这里的“高阶嵌入”分两种场景，对应不同实现方式：

• 场景1：将Riemannian度量网格嵌入到高维欧氏空间
  核心是让嵌入后的欧氏距离近似原度量下的测地距离：
  1. 基于顶点的Riemannian度量，计算顶点间的测地距离（可通过数值方法求解测地线，或用近似算法比如基于网格的最短路径，结合度量场调整边的权重）；
  2. 使用多维缩放（MDS）将测地距离矩阵映射到高维欧氏空间，CGAL本身没有MDS接口，但可以结合Eigen等线性代数库实现，再将嵌入后的坐标更新到网格顶点。
• 场景2：生成适配度量场的高阶网格（比如二次/三次三角剖分）
  1. 先用CGAL生成基础的度量适配三角剖分（如第一部分所述）；
  2. 使用CGAL的Subdivision_method_3模块进行高阶细分（比如Catmull-Clark或Loop细分）；
  3. 细分过程中，对新生成的顶点进行度量场插值（比如基于相邻顶点的度量矩阵做加权平均，权重参考细分权重），确保新顶点也带有适配的度量场。

关键注意事项 #

• CGAL的kernel是基础几何对象的实现，自定义度量不需要修改kernel，而是通过traits类扩展三角剖分的几何逻辑；
• Riemannian度量下的Delaunay剖分计算量远大于欧氏情况，若点集规模较大，需要优化测地圆判断的计算效率；
• 若不需要严格的Delaunay条件，也可以先生成欧氏Delaunay剖分，再基于度量场用CGAL的Mesh_optimization模块调整顶点位置，使网格适配度量场。

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