7.3 约束三角化 (Constrained_triangulation_2)

上一节: 7.2 二维Delaunay三角化下一节: 7.4 正则三角化

0. 直观理解：道路网络与地图绘制

0.1 生活类比：城市规划中的道路网络

想象你正在设计一座城市的道路网络。城市中有一些必须保留的原始道路（比如古老的商业街、高速公路），这些是城市的”约束”。无论城市如何发展，这些道路都不能被拆除或改变。

约束三角化就像是城市规划：

点：城市中的地标建筑、重要设施
约束边：必须保留的道路、铁路线
三角形：城市街区的划分

原始道路网络（约束边）        约束三角化结果
    
    A-------B                    A-------B
    |       |                    |\     /|
    |       道路1    -->         |  \ /  |
    |       |                    |  / \  |
    C-------D                    |/     \|
                                 C-------D
                                  道路2

0.2 为什么需要约束三角化？

场景1：地形建模

山脉的山脊线和山谷线必须保留
河流的走向不能改变
这些自然地物就是”约束边”

场景2：CAD设计

机械零件的边缘必须精确表示
不同材料的边界不能模糊
约束边确保设计精度

场景3：有限元分析

材料裂纹必须作为约束边
不同材料的分界面必须保留
否则会导致计算错误

0.3 约束Delaunay vs 严格Delaunay

严格Delaunay三角化              约束Delaunay三角化
（不允许约束边）                 （允许约束边）

    A-------B                      A-------B
     \     /                        \     /
      \   /                          \   /
       \ /                            \ /
        C                              C
        |                              |
        |                              |
        D                              D
        
Delaunay性质：                    约束边AB强制保留：
- 所有边都满足空圆性质            - AB边可能不满足Delaunay
- 但可能丢失重要特征              - 但保留了关键几何特征

关键区别：

严格Delaunay：所有边都满足空圆性质，但可能丢失几何特征
约束Delaunay：强制保留约束边，在约束存在的前提下尽可能满足Delaunay性质

1. 数学理论基础

1.1 约束三角化的定义

约束三角化 (Constrained Triangulation) 是在给定平面点集和一组约束边 (Constrained Edges) 的条件下构建的三角化，要求：

所有约束边必须作为三角化的边出现
三角化在其他方面保持标准性质（无重叠、覆盖凸包）

形式化定义：给定：

点集 $P = {p_{1}, p_{2}, ..., p_{n}}$
约束边集 $C = {(a_{1}, b_{1}), (a_{2}, b_{2}), ..., (a_{m}, b_{m})}$ ，其中 $a_{i}, b_{i} \in P$

约束三角化 $CT (P, C)$ 满足： $\forall (a, b) \in C : (a, b) \in edges (CT (P, C))$

1.2 约束边处理过程可视化

约束边插入的完整过程：

步骤1: 初始三角化（无约束）
        A
       /|\
      / | \
     /  |  \
    /   |   \
   /    |    \
  B-----C-----D

步骤2: 插入约束边AD（穿过现有边）
        A
       /|\
      / | \
     /  |  \   <-- 约束边AD（红色）
    /   |   \      与边BC相交
   /    |    \
  B-----C-----D

步骤3: 移除相交边，形成多边形洞
        A
       /|\
      / | \
     /  +  \   <-- 交点E
    / /   \ \
   / /     \ \
  B---------D
      洞区域

步骤4: 重新三角化洞区域
        A
       /|\
      / | \
     E--|--E    <-- 新顶点E
    / \ | / \
   /   \|/   \
  B-----C-----D

步骤5: 标记约束边
        A
       /|\
      / | \
     #==|==#    <-- 约束边（标记为#）
    / \ | / \
   /   \|/   \
  B-----C-----D

1.3 约束边与约束Delaunay

1.3.1 纯约束三角化

仅要求约束边存在，不优化三角形质量。可能产生细长三角形。

1.3.2 约束Delaunay三角化 (Constrained Delaunay Triangulation, CDT)

在满足约束的前提下，尽可能满足Delaunay性质：

一个三角形是约束Delaunay的，如果其外接圆内部不包含任何可见的其他点（约束边阻挡视线）。

形式化定义：三角形 $t$ 是约束Delaunay的，如果对于外接圆内的任意其他点 $p$ ，线段 $pq$ 与某条约束边相交（ $q$ 是 $t$ 的任意顶点）。

可见性概念可视化：

约束边阻挡视线示例：

    p (点在外接圆内)
     \
      \   约束边
       \   |
        \  |  <-- 约束边阻挡了p与三角形的视线
         \ |
          \|
    A------B------C
          三角形
          
结果：虽然p在外接圆内，但由于约束边阻挡，
      三角形ABC仍然是约束Delaunay的

1.4 约束边相交处理

当约束边相交时，需要处理交点：

相交类型	处理方式
在端点相交	允许，共享顶点
在内部相交	需要分割约束边
重叠	合并为一条边

相交处理可视化：

情况1: 端点相交（允许）
    A-----B
          |
          C-----D
    
情况2: 内部相交（需要分割）
    A-----+-----B  -->  A-----E-----B
          |E                |
    C-----+-----D           E
                            |
                      C-----D
    
情况3: 重叠（合并）
    A-----B-----C  -->  A-----------C
          B

CGAL提供三种策略（通过Itag模板参数）：

// 1. No_intersection_tag - 假设约束边不相交
// 2. Exact_intersections_tag - 使用精确计算处理交点
// 3. Exact_predicates_tag - 使用精确谓词，近似构造

1.5 多边形区域三角化

约束三角化的重要应用是将多边形区域（可能含洞）分解为三角形：

多边形边界作为约束边 -> 约束三角化 -> 删除外部三角形 -> 多边形三角化

含洞多边形处理过程：

步骤1: 定义外边界（逆时针）
        A-----------B
        |           |
        |    洞     |
        |           |
        D-----------C

步骤2: 定义洞边界（顺时针）
        A-----------B
        |   E-----F |
        |   | 洞  | |
        |   H-----G |
        D-----------C

步骤3: 约束三角化
        A-----------B
        | \   1   / |
        |   E---F   |
        | 2 |\3/| 4 |
        |   H---G   |
        | /   5   \ |
        D-----------C

步骤4: 标记内部/外部
        A-----------B
        | \   1   / |
        |   E---F   |
        | 2 |\3/| 4 |   <-- 1,2,3,4,5为内部三角形
        |   H---G   |
        | /   5   \ |
        D-----------C

2. 类架构分析

2.1 继承层次

Triangulation_2<Traits, Tds>
    └── Constrained_triangulation_2<Traits, Tds, Itag>
            └── Constrained_Delaunay_triangulation_2<Traits, Tds, Itag>
            └── Constrained_triangulation_plus_2<Traits, Tds, Itag>  // 支持边移除

2.2 模板参数

template <class Traits_, class Tds_ = Default, class Itag_ = Default>
class Constrained_triangulation_2 : public Triangulation_2<Traits_, Tds_> {
    // ...
};

Itag (Intersection Tag)

// 1. No_intersection_tag
//    - 用户保证约束边不相交（除端点外）
//    - 最高效，无运行时开销
 
// 2. Exact_predicates_tag  
//    - 使用精确谓词检测相交
//    - 交点用双精度近似表示
//    - 适用于大多数应用场景
 
// 3. Exact_intersections_tag
//    - 使用精确数表示交点
//    - 无精度损失，但速度较慢
//    - 适用于需要精确结果的场景

2.3 核心类型与方法

// 继承自Triangulation_2的所有类型
using Base = Triangulation_2<Traits, Tds>;
 
// 约束相关类型
typedef std::pair<Vertex_handle, Vertex_handle> Constraint;
typedef std::list<Edge> Constraint_list;
 
// 核心方法
// 插入约束边
void insert_constraint(Point a, Point b);  // 通过点
void insert_constraint(Vertex_handle va, Vertex_handle vb);  // 通过顶点句柄
void insert_constraint(Iterator first, Iterator last);  // 批量插入
 
// 约束边查询
bool is_constrained(Edge e);           // 边是否为约束边
bool is_constrained(Vertex_handle v);  // 顶点是否关联约束边
int number_of_constraints();           // 约束边数量
 
// 约束边遍历
Constraint_iterator constraints_begin();
Constraint_iterator constraints_end();
 
// 与子约束相关的操作 (Constrained_triangulation_plus_2)
void remove_constraint(Vertex_handle va, Vertex_handle vb);  // 移除约束
void split_constraint(Vertex_handle va, Vertex_handle vb, Vertex_handle vc);  // 分割

2.4 内部表示

// 面基类扩展，添加约束边标记
template <class Gt, class Fb = Triangulation_face_base_2<Gt>>
class Constrained_triangulation_face_base_2 : public Fb {
private:
    bool constraint_flags[3];  // 标记三条边是否为约束边
    // 边i连接vertices[(i+1)%3]和vertices[(i+2)%3]
    
public:
    bool is_constrained(int i) const { return constraint_flags[i]; }
    void set_constraint(int i, bool b) { constraint_flags[i] = b; }
    void set_constraints(bool c0, bool c1, bool c2) {
        constraint_flags[0] = c0;
        constraint_flags[1] = c1;
        constraint_flags[2] = c2;
    }
};

3. 实现细节

3.1 约束边插入算法

// 插入约束边的主算法
void insert_constraint(Vertex_handle va, Vertex_handle vb) {
    // 1. 检查边是否已存在
    if (are_there_incident_constraints(va, vb)) {
        return;  // 约束已存在
    }
    
    // 2. 找到与线段va-vb相交的所有边
    std::list<Edge> intersected_edges;
    find_intersected_edges(va, vb, intersected_edges);
    
    // 3. 处理相交边（根据Itag策略）
    process_intersections(va, vb, intersected_edges);
    
    // 4. 移除相交边，创建"多边形洞"
    std::list<Face_handle> hole;
    make_hole(va, vb, intersected_edges, hole);
    
    // 5. 重新三角化洞区域
    triangulate_hole(hole);
    
    // 6. 标记约束边
    mark_constraints(va, vb);
}
 
// 查找与线段相交的所有边
void find_intersected_edges(Vertex_handle va, Vertex_handle vb,
                            std::list<Edge>& edges) {
    // 使用行走法找到所有相交边
    Face_handle start = locate(va->point());
    
    // 从va向vb遍历，收集相交边
    walk_and_collect(va->point(), vb->point(), start, edges);
}
 
// 创建多边形洞
void make_hole(Vertex_handle va, Vertex_handle vb,
               const std::list<Edge>& intersected_edges,
               std::list<Face_handle>& hole) {
    // 删除所有与约束边相交的边
    // 形成一个多边形空洞
    for (const Edge& e : intersected_edges) {
        remove_edge(e);
    }
}
 
// 重新三角化洞
void triangulate_hole(std::list<Face_handle>& hole) {
    // 使用耳切割法或递归分割
    // 将多边形分解为三角形
}

3.2 约束Delaunay算法

// Constrained_Delaunay_triangulation_2的边翻转
void restore_constrained_delaunay() {
    std::stack<Edge> edges_to_check;
    
    // 收集所有内部边
    for (auto eit = edges_begin(); eit \!= edges_end(); ++eit) {
        if (\!is_constrained(*eit) && \!is_infinite(*eit)) {
            edges_to_check.push(*eit);
        }
    }
    
    // 边翻转循环
    while (\!edges_to_check.empty()) {
        Edge e = edges_to_check.top();
        edges_to_check.pop();
        
        // 跳过约束边
        if (is_constrained(e)) continue;
        
        Face_handle f = e.first;
        int i = e.second;
        Face_handle n = f->neighbor(i);
        
        if (is_infinite(n)) continue;
        
        // 检查是否可翻转且需要翻转
        if (is_flipable(f, i) && needs_flip(f, i)) {
            // 执行翻转
            flip(f, i);
            
            // 将新边加入检查队列
            edges_to_check.push(Edge(f, (i+1)%3));
            edges_to_check.push(Edge(n, (n->index(f)+1)%3));
        }
    }
}
 
// 检查是否需要翻转（考虑约束）
bool needs_flip(Face_handle f, int i) {
    Face_handle n = f->neighbor(i);
    Vertex_handle v_op = n->vertex(n->index(f));
    
    // 标准Delaunay测试
    if (side_of_oriented_circle(f, v_op->point()) \!= ON_POSITIVE_SIDE) {
        return false;  // 已经满足Delaunay
    }
    
    // 检查翻转后的边是否会与约束边相交
    // 如果会相交，则不能翻转
    Vertex_handle v1 = f->vertex((i+1)%3);
    Vertex_handle v2 = f->vertex((i+2)%3);
    
    // 新边将是v_op到f中对面边i的顶点
    // 需要检查这条新边是否与任何约束边相交
    
    return can_flip_without_violating_constraints(f, i);
}

3.3 交点处理 (Exact_intersections_tag)

// 当约束边相交时，需要插入交点
template <class Gt>
class Exact_intersections_tag {
    // 使用精确数表示交点坐标
    using FT = Gt::FT;  // 可能是Gmpq等精确数类型
    
    Point compute_intersection(const Segment& s1, const Segment& s2) {
        // 使用精确计算求交点
        // 返回精确坐标
    }
};
 
// 插入带交点的约束
void insert_constraint_with_intersections(Point a, Point b) {
    // 1. 找到所有交点
    std::vector<Point> intersections;
    find_all_intersections(a, b, intersections);
    
    // 2. 按在线段上的顺序排序
    std::sort(intersections.begin(), intersections.end(), 
              [&](const Point& p1, const Point& p2) {
                  return squared_distance(a, p1) < squared_distance(a, p2);
              });
    
    // 3. 逐段插入约束
    Point prev = a;
    for (const Point& inter : intersections) {
        insert_constraint(prev, inter);
        prev = inter;
    }
    insert_constraint(prev, b);
}

4. 渐进式示例教程

4.1 最简示例：单条约束边

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <iostream>
 
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using CDT = CGAL::Constrained_Delaunay_triangulation_2<K>;
using Point = K::Point_2;
 
int main() {
    CDT cdt;
    
    // 插入4个点形成矩形
    Point a(0, 0), b(10, 0), c(10, 10), d(0, 10);
    cdt.insert(a);
    cdt.insert(b);
    cdt.insert(c);
    cdt.insert(d);
    
    // 插入一条对角线约束
    cdt.insert_constraint(a, c);
    
    std::cout << "顶点数: " << cdt.number_of_vertices() << std::endl;
    std::cout << "约束边数: " << cdt.number_of_constraints() << std::endl;
    
    return 0;
}

4.2 基础约束三角化

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_triangulation_2.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <iostream>
#include <vector>
 
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
 
// 使用Exact_predicates_tag处理可能的相交
using CDT = CGAL::Constrained_Delaunay_triangulation_2<
    K,
    CGAL::Default,
    CGAL::Exact_predicates_tag
>;
 
using Point = K::Point_2;
using Vertex_handle = CDT::Vertex_handle;
using Edge = CDT::Edge;
 
int main() {
    CDT cdt;
    
    // ========== 插入点 ==========
    std::cout << "=== 约束Delaunay三角化 ===" << std::endl;
    
    // 创建一个矩形区域的点
    std::vector<Point> points = {
        Point(0, 0), Point(10, 0), Point(10, 10), Point(0, 10),  // 角点
        Point(5, 0), Point(10, 5), Point(5, 10), Point(0, 5),    // 边中点
        Point(5, 5)                                              // 中心
    };
    
    // 批量插入点
    std::vector<Vertex_handle> vertices;
    for (const auto& p : points) {
        vertices.push_back(cdt.insert(p));
    }
    
    std::cout << "插入 " << cdt.number_of_vertices() << " 个顶点" << std::endl;
    
    // ========== 插入约束边 ==========
    std::cout << "\n=== 插入约束边 ===" << std::endl;
    
    // 约束1: 矩形外边界
    cdt.insert_constraint(vertices[0], vertices[4]);  // (0,0) -> (5,0)
    cdt.insert_constraint(vertices[4], vertices[1]);  // (5,0) -> (10,0)
    cdt.insert_constraint(vertices[1], vertices[5]);  // (10,0) -> (10,5)
    cdt.insert_constraint(vertices[5], vertices[2]);  // (10,5) -> (10,10)
    cdt.insert_constraint(vertices[2], vertices[6]);  // (10,10) -> (5,10)
    cdt.insert_constraint(vertices[6], vertices[3]);  // (5,10) -> (0,10)
    cdt.insert_constraint(vertices[3], vertices[7]);  // (0,10) -> (0,5)
    cdt.insert_constraint(vertices[7], vertices[0]);  // (0,5) -> (0,0)
    
    std::cout << "外边界约束插入完成" << std::endl;
    
    // 约束2: 对角线（将矩形分成两个三角形区域）
    cdt.insert_constraint(vertices[0], vertices[2]);  // (0,0) -> (10,10)
    
    std::cout << "对角线约束插入完成" << std::endl;
    
    // 约束3: 内部线段
    cdt.insert_constraint(vertices[4], vertices[6]);  // (5,0) -> (5,10) - 垂直线
    
    std::cout << "内部约束插入完成" << std::endl;
    
    // ========== 统计信息 ==========
    std::cout << "\n=== 统计信息 ===" << std::endl;
    std::cout << "顶点数: " << cdt.number_of_vertices() << std::endl;
    std::cout << "面数: " << cdt.number_of_faces() << std::endl;
    std::cout << "约束边数: " << cdt.number_of_constraints() << std::endl;
    
    // ========== 遍历约束边 ==========
    std::cout << "\n=== 约束边列表 ===" << std::endl;
    int constraint_count = 0;
    for (auto eit = cdt.edges_begin(); eit \!= cdt.edges_end(); ++eit) {
        if (cdt.is_constrained(*eit)) {
            auto seg = cdt.segment(*eit);
            std::cout << "约束边 " << ++constraint_count << ": ("
                      << seg.source().x() << ", " << seg.source().y() << ") -> ("
                      << seg.target().x() << ", " << seg.target().y() << ")" << std::endl;
        }
    }
    
    // ========== 验证约束 ==========
    std::cout << "\n=== 约束验证 ===" << std::endl;
    
    // 检查特定边是否为约束边
    for (auto eit = cdt.finite_edges_begin(); eit \!= cdt.finite_edges_end(); ++eit) {
        Edge e = *eit;
        if (cdt.is_constrained(e)) {
            // 验证这条边确实存在于三角化中
            std::cout << "约束边验证通过" << std::endl;
            break;
        }
    }
    
    // ========== 遍历所有三角形 ==========
    std::cout << "\n=== 三角形列表 ===" << std::endl;
    int triangle_count = 0;
    for (auto fit = cdt.finite_faces_begin(); fit \!= cdt.finite_faces_end(); ++fit) {
        std::cout << "三角形 " << ++triangle_count << ": ";
        for (int i = 0; i < 3; ++i) {
            Point p = fit->vertex(i)->point();
            std::cout << "(" << p.x() << ", " << p.y() << ") ";
        }
        
        // 标记约束边
        std::cout << "[约束边: ";
        for (int i = 0; i < 3; ++i) {
            if (fit->is_constrained(i)) {
                std::cout << i << " ";
            }
        }
        std::cout << "]" << std::endl;
    }
    
    return 0;
}

4.3 多边形三角化

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/mark_domain_in_triangulation.h>
#include <iostream>
#include <vector>
 
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using CDT = CGAL::Constrained_Delaunay_triangulation_2<K>;
 
using Point = K::Point_2;
using Vertex_handle = CDT::Vertex_handle;
using Face_handle = CDT::Face_handle;
 
// 多边形三角化类
class PolygonTriangulation {
public:
    // 添加多边形外边界（逆时针）
    void add_outer_boundary(const std::vector<Point>& polygon) {
        if (polygon.size() < 3) return;
        
        // 插入点
        std::vector<Vertex_handle> vertices;
        for (const auto& p : polygon) {
            vertices.push_back(cdt.insert(p));
        }
        
        // 插入约束边
        for (size_t i = 0; i < vertices.size(); ++i) {
            size_t j = (i + 1) % vertices.size();
            cdt.insert_constraint(vertices[i], vertices[j]);
        }
    }
    
    // 添加洞（顺时针）
    void add_hole(const std::vector<Point>& hole) {
        if (hole.size() < 3) return;
        
        std::vector<Vertex_handle> vertices;
        for (const auto& p : hole) {
            vertices.push_back(cdt.insert(p));
        }
        
        for (size_t i = 0; i < vertices.size(); ++i) {
            size_t j = (i + 1) % vertices.size();
            cdt.insert_constraint(vertices[i], vertices[j]);
        }
    }
    
    // 获取多边形内部的三角形
    std::vector<Face_handle> get_interior_triangles() {
        // 标记域（区分内部和外部）
        mark_domains();
        
        std::vector<Face_handle> interior;
        for (auto fit = cdt.finite_faces_begin(); fit \!= cdt.finite_faces_end(); ++fit) {
            if (fit->info().in_domain) {  // 需要自定义面信息
                interior.push_back(fit);
            }
        }
        return interior;
    }
    
    CDT& get_cdt() { return cdt; }
    
private:
    CDT cdt;
    
    // 使用BFS标记域
    void mark_domains() {
        // 默认所有面为外部
        for (auto fit = cdt.finite_faces_begin(); fit \!= cdt.finite_faces_end(); ++fit) {
            fit->info().in_domain = false;
        }
        
        // 从无限面开始BFS，标记可达面为外部
        std::queue<Face_handle> queue;
        std::set<Face_handle> visited;
        
        // 找到无限面
        Face_handle infinite_face = cdt.infinite_face();
        queue.push(infinite_face);
        visited.insert(infinite_face);
        
        while (\!queue.empty()) {
            Face_handle f = queue.front();
            queue.pop();
            
            for (int i = 0; i < 3; ++i) {
                Face_handle n = f->neighbor(i);
                
                // 如果邻接面未访问且不是约束边，则可达
                if (visited.find(n) == visited.end() && \!cdt.is_constrained(Edge(f, i))) {
                    visited.insert(n);
                    queue.push(n);
                }
            }
        }
        
        // 未访问的有限面为内部
        for (auto fit = cdt.finite_faces_begin(); fit \!= cdt.finite_faces_end(); ++fit) {
            if (visited.find(fit) == visited.end()) {
                fit->info().in_domain = true;
            }
        }
    }
};
 
// 自定义面信息
struct FaceInfo {
    bool in_domain;
    FaceInfo() : in_domain(false) {}
};
 
// 重新定义带信息的CDT
using Fb = CGAL::Triangulation_face_base_with_info_2<FaceInfo, K>;
using Tds = CGAL::Triangulation_data_structure_2<
    CGAL::Triangulation_vertex_base_2<K>,
    Fb
>;
using CDT_with_info = CGAL::Constrained_Delaunay_triangulation_2<K, Tds>;
 
int main() {
    CDT_with_info cdt;
    
    std::cout << "=== 多边形三角化 ===" << std::endl;
    
    // 定义外边界（逆时针）- 矩形
    std::vector<Point> outer = {
        Point(0, 0), Point(10, 0), Point(10, 10), Point(0, 10)
    };
    
    // 定义洞（顺时针）- 内部矩形
    std::vector<Point> hole = {
        Point(3, 3), Point(3, 7), Point(7, 7), Point(7, 3)
    };
    
    // 插入外边界点
    std::vector<Vertex_handle> outer_vertices;
    for (const auto& p : outer) {
        outer_vertices.push_back(cdt.insert(p));
    }
    
    // 插入外边界约束
    for (size_t i = 0; i < outer_vertices.size(); ++i) {
        size_t j = (i + 1) % outer_vertices.size();
        cdt.insert_constraint(outer_vertices[i], outer_vertices[j]);
    }
    
    std::cout << "外边界插入完成" << std::endl;
    
    // 插入洞点
    std::vector<Vertex_handle> hole_vertices;
    for (const auto& p : hole) {
        hole_vertices.push_back(cdt.insert(p));
    }
    
    // 插入洞约束
    for (size_t i = 0; i < hole_vertices.size(); ++i) {
        size_t j = (i + 1) % hole_vertices.size();
        cdt.insert_constraint(hole_vertices[i], hole_vertices[j]);
    }
    
    std::cout << "洞插入完成" << std::endl;
    
    // 标记域
    std::queue<Face_handle> queue;
    std::set<Face_handle> visited;
    
    // 从无限面开始BFS
    queue.push(cdt.infinite_face());
    visited.insert(cdt.infinite_face());
    
    while (\!queue.empty()) {
        Face_handle f = queue.front();
        queue.pop();
        
        for (int i = 0; i < 3; ++i) {
            Face_handle n = f->neighbor(i);
            if (visited.find(n) == visited.end()) {
                // 检查是否为约束边
                bool is_constrained = false;
                if (\!cdt.is_infinite(f)) {
                    is_constrained = f->is_constrained(i);
                }
                
                if (\!is_constrained) {
                    visited.insert(n);
                    queue.push(n);
                }
            }
        }
    }
    
    // 标记内部面
    for (auto fit = cdt.finite_faces_begin(); fit \!= cdt.finite_faces_end(); ++fit) {
        fit->info().in_domain = (visited.find(fit) == visited.end());
    }
    
    // 统计内部三角形
    int interior_count = 0;
    int exterior_count = 0;
    
    for (auto fit = cdt.finite_faces_begin(); fit \!= cdt.finite_faces_end(); ++fit) {
        if (fit->info().in_domain) {
            ++interior_count;
        } else {
            ++exterior_count;
        }
    }
    
    std::cout << "\n=== 三角化结果 ===" << std::endl;
    std::cout << "总面数: " << cdt.number_of_faces() << std::endl;
    std::cout << "内部三角形: " << interior_count << std::endl;
    std::cout << "外部三角形: " << exterior_count << std::endl;
    
    // 输出内部三角形
    std::cout << "\n=== 内部三角形 ===" << std::endl;
    int count = 0;
    for (auto fit = cdt.finite_faces_begin(); fit \!= cdt.finite_faces_end(); ++fit) {
        if (fit->info().in_domain) {
            std::cout << "三角形 " << ++count << ": ";
            for (int i = 0; i < 3; ++i) {
                Point p = fit->vertex(i)->point();
                std::cout << "(" << p.x() << ", " << p.y() << ") ";
            }
            std::cout << std::endl;
        }
    }
    
    return 0;
}

4.4 带交点处理的约束插入

#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <iostream>
#include <vector>
 
// 使用精确构造内核处理交点
using K = CGAL::Exact_predicates_exact_constructions_kernel;
using CDT = CGAL::Constrained_Delaunay_triangulation_2<
    K,
    CGAL::Default,
    CGAL::Exact_intersections_tag
>;
 
using Point = K::Point_2;
using Segment = K::Segment_2;
using Vertex_handle = CDT::Vertex_handle;
 
int main() {
    CDT cdt;
    
    std::cout << "=== 带交点处理的约束三角化 ===" << std::endl;
    
    // 插入初始点
    std::vector<Point> points = {
        Point(0, 0), Point(10, 0), Point(10, 10), Point(0, 10)
    };
    
    std::vector<Vertex_handle> vertices;
    for (const auto& p : points) {
        vertices.push_back(cdt.insert(p));
    }
    
    // 插入外边界
    for (size_t i = 0; i < vertices.size(); ++i) {
        size_t j = (i + 1) % vertices.size();
        cdt.insert_constraint(vertices[i], vertices[j]);
    }
    
    std::cout << "初始顶点数: " << cdt.number_of_vertices() << std::endl;
    
    // 插入两条相交的约束边
    // 约束1: 从(0,0)到(10,10)
    // 约束2: 从(0,10)到(10,0)
    // 这两条线在(5,5)相交
    
    std::cout << "\n插入第一条约束: (0,0) -> (10,10)" << std::endl;
    cdt.insert_constraint(points[0], points[2]);
    std::cout << "当前顶点数: " << cdt.number_of_vertices() << std::endl;
    
    std::cout << "\n插入第二条约束: (0,10) -> (10,0)" << std::endl;
    cdt.insert_constraint(points[3], points[1]);
    std::cout << "当前顶点数: " << cdt.number_of_vertices() << std::endl;
    
    // 由于Exact_intersections_tag，交点(5,5)被自动插入
    std::cout << "\n交点处理后顶点数: " << cdt.number_of_vertices() << std::endl;
    
    // 验证交点存在
    bool found_intersection = false;
    Point expected_intersection(5, 5);
    
    for (auto vit = cdt.finite_vertices_begin(); vit \!= cdt.finite_vertices_end(); ++vit) {
        // 检查是否有顶点在(5,5)
        if (vit->point() == expected_intersection) {
            found_intersection = true;
            std::cout << "找到交点: (5, 5)" << std::endl;
            break;
        }
    }
    
    if (\!found_intersection) {
        std::cout << "警告: 未找到预期交点" << std::endl;
    }
    
    // 输出所有顶点
    std::cout << "\n=== 所有顶点 ===" << std::endl;
    for (auto vit = cdt.finite_vertices_begin(); vit \!= cdt.finite_vertices_end(); ++vit) {
        Point p = vit->point();
        // 转换为double输出（精确数不能直接输出）
        double x = CGAL::to_double(p.x());
        double y = CGAL::to_double(p.y());
        std::cout << "(" << x << ", " << y << ")" << std::endl;
    }
    
    // 输出约束边
    std::cout << "\n=== 约束边 ===" << std::endl;
    for (auto eit = cdt.edges_begin(); eit \!= cdt.edges_end(); ++eit) {
        if (cdt.is_constrained(*eit)) {
            Segment seg = cdt.segment(*eit);
            double x1 = CGAL::to_double(seg.source().x());
            double y1 = CGAL::to_double(seg.source().y());
            double x2 = CGAL::to_double(seg.target().x());
            double y2 = CGAL::to_double(seg.target().y());
            std::cout << "(" << x1 << ", " << y1 << ") -> (" 
                      << x2 << ", " << y2 << ")" << std::endl;
        }
    }
    
    return 0;
}

4.5 实际应用：地形建模

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <fstream>
#include <iostream>
 
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using CDT = CGAL::Constrained_Delaunay_triangulation_2<K>;
using Point = K::Point_2;
 
// 地形特征线结构
struct TerrainFeature {
    std::string name;
    std::vector<Point> points;
    bool is_ridge;  // true=山脊, false=山谷
};
 
int main() {
    CDT cdt;
    
    std::cout << "=== 地形建模示例 ===" << std::endl;
    
    // 定义地形特征线
    std::vector<TerrainFeature> features = {
        // 山脊线
        {"Ridge1", {{10,10}, {30,25}, {50,30}, {70,25}, {90,10}}, true},
        // 山谷线
        {"Valley1", {{10,90}, {30,75}, {50,70}, {70,75}, {90,90}}, false},
        // 河流
        {"River", {{0,50}, {25,45}, {50,50}, {75,55}, {100,50}}, false}
    };
    
    // 插入地形点
    for (const auto& feature : features) {
        std::cout << "处理特征线: " << feature.name << std::endl;
        
        // 插入点并创建约束
        std::vector<CDT::Vertex_handle> vertices;
        for (const auto& p : feature.points) {
            vertices.push_back(cdt.insert(p));
        }
        
        // 创建约束边（确保特征线被保留）
        for (size_t i = 0; i < vertices.size() - 1; ++i) {
            cdt.insert_constraint(vertices[i], vertices[i+1]);
        }
    }
    
    // 插入随机地形点
    std::vector<Point> terrain_points = {
        {20, 20}, {40, 20}, {60, 20}, {80, 20},
        {20, 40}, {40, 40}, {60, 40}, {80, 40},
        {20, 60}, {40, 60}, {60, 60}, {80, 60},
        {20, 80}, {40, 80}, {60, 80}, {80, 80}
    };
    
    for (const auto& p : terrain_points) {
        cdt.insert(p);
    }
    
    std::cout << "\n地形三角化完成:" << std::endl;
    std::cout << "顶点数: " << cdt.number_of_vertices() << std::endl;
    std::cout << "三角形数: " << cdt.number_of_faces() << std::endl;
    std::cout << "约束边数: " << cdt.number_of_constraints() << std::endl;
    
    // 导出为VTK格式（用于可视化）
    std::ofstream vtk_file("terrain.vtk");
    vtk_file << "# vtk DataFile Version 3.0\n";
    vtk_file << "Terrain Triangulation\n";
    vtk_file << "ASCII\n";
    vtk_file << "DATASET UNSTRUCTURED_GRID\n";
    
    // 写入点
    vtk_file << "POINTS " << cdt.number_of_vertices() << " float\n";
    std::map<CDT::Vertex_handle, int> vertex_indices;
    int idx = 0;
    for (auto vit = cdt.finite_vertices_begin(); vit \!= cdt.finite_vertices_end(); ++vit) {
        Point p = vit->point();
        vtk_file << p.x() << " " << p.y() << " 0\n";
        vertex_indices[vit] = idx++;
    }
    
    // 写入三角形
    int num_triangles = 0;
    for (auto fit = cdt.finite_faces_begin(); fit \!= cdt.finite_faces_end(); ++fit) {
        ++num_triangles;
    }
    
    vtk_file << "CELLS " << num_triangles << " " << num_triangles * 4 << "\n";
    for (auto fit = cdt.finite_faces_begin(); fit \!= cdt.finite_faces_end(); ++fit) {
        vtk_file << "3 ";
        for (int i = 0; i < 3; ++i) {
            vtk_file << vertex_indices[fit->vertex(i)] << " ";
        }
        vtk_file << "\n";
    }
    
    vtk_file << "CELL_TYPES " << num_triangles << "\n";
    for (int i = 0; i < num_triangles; ++i) {
        vtk_file << "5\n";  // VTK_TRIANGLE
    }
    
    vtk_file.close();
    std::cout << "地形数据已导出到 terrain.vtk" << std::endl;
    
    return 0;
}

5. 复杂度分析

5.1 时间复杂度

操作	平均复杂度	最坏复杂度	说明
insert(p)	O(log n)	O(n)	点插入
insert_constraint(a,b)	O(n)	O(n²)	约束插入，取决于相交边数
批量约束插入	O(n²)	O(n²)	大量相交约束
is_constrained(e)	O(1)	O(1)	常数时间查询
遍历约束边	O(n)	O(n)	线性遍历
多边形三角化	O(n log n)	O(n²)	含洞的多边形

5.2 空间复杂度

存储项目	空间
基础三角化	O(n)
约束标记	O(n)
交点（Exact_intersections_tag）	O(k)，k为交点数
总计	O(n + k)

5.3 Itag选择指南

Itag	适用场景	性能	精度
No_intersection_tag	约束边不相交	最快	无交点
Exact_predicates_tag	一般应用	中等	近似交点
Exact_intersections_tag	需要精确结果	较慢	精确交点

6. 应用场景与参数选择

6.1 地理信息系统 (GIS)

// 地形建模中的特征线保留
// 山脊线、山谷线作为约束边

参数建议：

使用 Exact_predicates_tag（平衡性能和精度）
特征线采样密度：每10-50米一个点
约束边相交处理：自动（允许道路交叉口）

6.2 有限元分析

// 网格生成
// 材料边界、裂纹作为约束边

参数建议：

使用 Exact_intersections_tag（确保几何精确性）
材料边界必须严格保留
裂纹尖端需要特殊处理（加密）

6.3 计算机图形学

// 多边形剖分用于渲染
// 保留纹理边界

参数建议：

使用 No_intersection_tag（性能优先）
预处理确保约束边不相交
使用约束Delaunay优化三角形质量

6.4 路径规划

// 障碍物边界作为约束
// 在自由空间内规划路径

参数建议：

障碍物边界作为约束边
使用CDT构建导航网格
在三角形中心或边上进行路径规划

7. 总结

Constrained_triangulation_2及其Delaunay变体提供了：

约束边保证：强制特定边存在于三角化中
灵活的相交处理：三种策略适应不同需求
多边形处理能力：支持含洞的多边形三角化
质量优化：约束Delaunay尽可能保持三角形质量

约束三角化是连接纯几何三角化与实际工程应用的关键工具，广泛应用于GIS、CAD、有限元等领域。

Chunibyo

Explorer

7.3 约束三角化 (Constrained_triangulation_2)

7.3 约束三角化 (Constrained_triangulation_2) §

0. 直观理解：道路网络与地图绘制 §

0.1 生活类比：城市规划中的道路网络 §

0.2 为什么需要约束三角化？ §

0.3 约束Delaunay vs 严格Delaunay §

1. 数学理论基础 §

1.1 约束三角化的定义 §

1.2 约束边处理过程可视化 §

1.3 约束边与约束Delaunay §

1.3.1 纯约束三角化 §

1.3.2 约束Delaunay三角化 (Constrained Delaunay Triangulation, CDT) §

1.4 约束边相交处理 §

1.5 多边形区域三角化 §

2. 类架构分析 §

2.1 继承层次 §

2.2 模板参数 §

Itag (Intersection Tag) §

2.3 核心类型与方法 §

2.4 内部表示 §

3. 实现细节 §

3.1 约束边插入算法 §

3.2 约束Delaunay算法 §

3.3 交点处理 (Exact_intersections_tag) §

4. 渐进式示例教程 §

4.1 最简示例：单条约束边 §

4.2 基础约束三角化 §

4.3 多边形三角化 §

4.4 带交点处理的约束插入 §

4.5 实际应用：地形建模 §

5. 复杂度分析 §

5.1 时间复杂度 §

5.2 空间复杂度 §

5.3 Itag选择指南 §

6. 应用场景与参数选择 §

6.1 地理信息系统 (GIS) §

6.2 有限元分析 §

6.3 计算机图形学 §

6.4 路径规划 §

7. 总结 §

Graph View

Table of Contents

Backlinks

7.3 约束三角化 (Constrained_triangulation_2)

0. 直观理解：道路网络与地图绘制

0.1 生活类比：城市规划中的道路网络

0.2 为什么需要约束三角化？

0.3 约束Delaunay vs 严格Delaunay

1. 数学理论基础

1.1 约束三角化的定义

1.2 约束边处理过程可视化

1.3 约束边与约束Delaunay

1.3.1 纯约束三角化

1.3.2 约束Delaunay三角化 (Constrained Delaunay Triangulation, CDT)

1.4 约束边相交处理

1.5 多边形区域三角化

2. 类架构分析

2.1 继承层次

2.2 模板参数

Itag (Intersection Tag)

2.3 核心类型与方法

2.4 内部表示

3. 实现细节

3.1 约束边插入算法

3.2 约束Delaunay算法

3.3 交点处理 (Exact_intersections_tag)

4. 渐进式示例教程

4.1 最简示例：单条约束边

4.2 基础约束三角化

4.3 多边形三角化

4.4 带交点处理的约束插入

4.5 实际应用：地形建模

5. 复杂度分析

5.1 时间复杂度

5.2 空间复杂度

5.3 Itag选择指南

6. 应用场景与参数选择

6.1 地理信息系统 (GIS)

6.2 有限元分析

6.3 计算机图形学

6.4 路径规划

7. 总结