CGAL如何实现类似Demo的特征折线约束表面网格分割

现有包含多处尖锐特征的表面网格，需要以尖锐特征构成的折线为边界分割网格。CGAL Demo中的Detect Sharp Features功能可实现该效果，效果见附图。目前已能通过domain.detect_features()接口获取特征折线，该接口内部会调用add_features_from_split_graph_into_polylines()方法，但无法提取折线围合区域内的表面面片、完成网格切割。
CGAL Demo效果示例

现有代码 #

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Polyhedral_mesh_domain_with_features_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/IO/output_to_vtu.h>

#include <CGAL/IO/facets_in_complex_3_to_triangle_mesh.h>
#include <CGAL/Surface_mesh.h>


// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Mesh_polyhedron_3<K>::type Polyhedron;
typedef CGAL::Polyhedral_mesh_domain_with_features_3<K> Mesh_domain;
#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif
// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain, CGAL::Default, Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<
    Tr, Mesh_domain::Corner_index, Mesh_domain::Curve_index> C3t3;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
// To avoid verbose function and named parameters call
using namespace CGAL::parameters;

//Mesh
typedef CGAL::Surface_mesh<K::Point_3> Mesh;

int main(int argc, char*argv[])
{
    const char* fname = (argc > 1) ? argv[1] : "data/fandisk.off";
    std::ifstream input(fname);
    Polyhedron polyhedron;
    input >> polyhedron;
    if (input.fail()) {
        std::cerr << "Error: Cannot read file " << fname << std::endl;
        return EXIT_FAILURE;
    }
    if (!CGAL::is_triangle_mesh(polyhedron)) {
        std::cerr << "Input geometry is not triangulated." << std::endl;
        return EXIT_FAILURE;
    }
    // Create domain
    Mesh_domain domain(polyhedron);
    // Get sharp features
    domain.detect_features();
    // Mesh criteria
    Mesh_criteria criteria(edge_size = 0.025,
        facet_angle = 25, facet_size = 0.05, facet_distance = 0.005,
        cell_radius_edge_ratio = 3, cell_size = 0.05);
    // Mesh generation
    C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
    
    // use cgal demo open vtu file can't get desired result
    // Output
    /*std::ofstream file("out.vtu");
    CGAL::output_to_vtu(file, c3t3, CGAL::IO::ASCII);*/

    //use function below can get mesh but doesn't contain segmentation info
    Mesh mesh;
    facets_in_complex_3_to_triangle_mesh(c3t3, mesh);
    std::ofstream output("mesh_smoothed.off");
    output.precision(17);
    output << mesh;


    // Could be replaced by:
    // c3t3.output_to_medit(file);
    return EXIT_SUCCESS;
}

当前代码的问题是用了facets_in_complex_3_to_triangle_mesh()直接导出网格，这个接口只会导出几何面片，不会保留特征边界分割后的面片分块标记。detect_features()执行后，生成网格时所有落在尖锐特征折线上的边都会被标记为约束边，表面网格会被这些边自然切分为多个连通区域，每个区域的面片都绑定了唯一的表面patch索引，直接读取这个索引就能拿到分块结果，不需要自己做区域围合计算。

核心实现逻辑 #

• 跳过facets_in_complex_3_to_triangle_mesh()调用，手动遍历C3t3结构中的表面面片构造输出网格
• 对每个表面面片，调用c3t3.surface_patch_index()获取所属的patch ID，同一个ID的面片就属于同一块被特征折线围合的区域
• 给输出网格添加面片属性存储patch ID，后续可以按ID分组提取子网格，或者直接把ID作为标量字段写入VTU文件，就能在CGAL Demo里看到分块效果

修改后的关键代码 #

替换原代码中facets_in_complex_3_to_triangle_mesh相关的导出逻辑即可：

Mesh mesh;
typedef boost::graph_traits<Mesh>::vertex_descriptor VertexDesc;
typedef Mesh::Property_map<Mesh::Face_index, int> FacePatchMap;

// 为输出网格添加面片分块ID属性
FacePatchMap patch_id_map = mesh.add_property_map<Mesh::Face_index, int>("f:patch_id", -1).first;
// 缓存三角剖分顶点到输出网格顶点的映射，避免重复插入顶点
std::map<Tr::Vertex_handle, VertexDesc> vertex_mapping;

// 遍历所有三角面片，仅保留属于表面的部分
for (auto fit = c3t3.facets_begin(); fit != c3t3.facets_end(); ++fit) {
    if (!c3t3.is_in_complex(*fit)) continue;

    typename Tr::Cell_handle cell = fit->first;
    int facet_on_cell_idx = fit->second;
    // 获取当前面片所属的分块ID
    const auto patch_idx = c3t3.surface_patch_index(cell, facet_on_cell_idx);
    const int block_id = static_cast<int>(patch_idx);

    // 构造当前三角面片的三个顶点
    std::vector<VertexDesc> tri_verts;
    for (int i = 0; i < 3; ++i) {
        const int v_local_idx = Tr::vertex_triple_index(facet_on_cell_idx, i);
        auto tr_vh = cell->vertex(v_local_idx);
        if (!vertex_mapping.count(tr_vh)) {
            vertex_mapping[tr_vh] = mesh.add_vertex(tr_vh->point());
        }
        tri_verts.push_back(vertex_mapping[tr_vh]);
    }

    // 向输出网格添加面片，写入分块ID
    auto new_face = mesh.add_face(tri_verts[0], tri_verts[1], tri_verts[2]);
    if (new_face != Mesh::null_face()) {
        patch_id_map[new_face] = block_id;
    }
}

// 导出带分块属性的网格，OFF格式会自动保留自定义属性
std::ofstream output("mesh_segmented.off");
output.precision(17);
output << mesh;

补充说明 #

• 如果不需要对原始网格做重网格化，只是想给输入网格做尖锐边分割，不需要走make_mesh_3的体网格生成流程，直接调用CGAL::Polygon_mesh_processing::detect_sharp_edges()提取尖锐边，再用带边约束的连通域计算接口即可，运行效率高很多。
• 如果需要把每个分块导出为独立网格，遍历所有面片按block_id分组，分别构造子网格即可。
• 之前导出的VTU文件看不到分块效果，是因为没有把patch ID作为单元数据写入VTU，手动添加标量字段即可正常显示。

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