Chunibyo

Search

SearchSearch

7.4 正则三角化 (Regular_triangulation_2)

Apr 04, 2026, 18 min read

7.4 正则三角化 (Regular_triangulation_2)

上一节: 7.2 二维Delaunay三角化 下一节: 7.5 特殊三角化

1. 数学理论基础

1.1 加权点与幂距离

加权点 (Weighted Point) 是带有权重的平面点：

其中 是位置， 是权重。

幂距离 (Power Distance) 定义为：

几何解释：

  • 当 ：加权点对应以 为中心、 为半径的圆
  • 当 ：退化为普通点
  • 当 ：对应虚圆

1.2 幂图 (Power Diagram)

幂图是Voronoi图在加权点上的推广：

幂图的性质：

  • 单元可能是空集（被隐藏的加权点）
  • 单元是凸多边形（可能无界）
  • 边的数目可能少于Voronoi图

1.3 正则三角化的定义

正则三角化 (Regular Triangulation) 是幂图的对偶图：

两个加权点的单元共享一条边，当且仅当它们在正则三角化中通过一条边连接。

形式化定义：对于加权点集 ，正则三角化 满足：

1.4 空外接圆性质的推广

正则三角化推广了Delaunay的空外接圆性质：

正交圆测试 (Orthogonality Test)：

三角形 是正则的，如果不存在其他加权点 使得：

\text{orthogonal_circle}(\hat{p}_i, \hat{p}_j, \hat{p}_k) \text{ 包含 } p_l \text{ 且 } w_l > \text{power}(p_l, \text{circle})

等价表述：对于三角形的外接正交圆 ，所有其他加权点的中心都在 外部，或者在 上但权重更小。

1.5 隐藏点 (Hidden Points)

当加权点的权重过大时，其幂单元可能为空，该点被隐藏：

即存在另一个加权点的圆完全包含 。

隐藏点的性质：

  • 隐藏点不参与正则三角化
  • 隐藏点可能被”揭示”（当其他点被删除时）
  • 正则三角化只包含可见点

1.6 与Delaunay的关系

特性DelaunayRegular
输入普通点加权点
对偶Voronoi图幂图
空圆测试外接圆正交圆
隐藏点无可能有
退化情况四点共圆四点共正交圆

当所有权重为0时，正则三角化退化为Delaunay三角化。

2. 类架构分析

2.1 继承层次

Triangulation_2<Traits, Tds>
    └── Delaunay_triangulation_2<Traits, Tds>
            └── Regular_triangulation_2<Traits, Tds>

注意：Regular_triangulation_2继承自Delaunay_triangulation_2，而非直接继承Triangulation_2。

2.2 模板参数

template <class Traits_, class Tds_ = Default>
class Regular_triangulation_2 : public Delaunay_triangulation_2<Traits_, Tds_> {
    // ...
};

Traits要求

class RegularTriangulationTraits_2 : public DelaunayTriangulationTraits_2 {
public:
    // 加权点类型
    typedef Weighted_point_2 Weighted_point;
    
    // 构造加权点
    typedef Construct_weighted_point_2 Construct_weighted_point;
    
    // 幂距离
    typedef Power_distance_2 Power_distance;
    
    // 正交圆测试
    typedef Power_test_2 Power_test;  // 相当于Delaunay的side_of_oriented_circle
    
    // 从加权点提取普通点
    typedef Construct_point_2 Construct_point;
};

2.3 核心类型与方法

// 继承自Delaunay_triangulation_2的所有类型
using Base = Delaunay_triangulation_2<Traits, Tds>;
 
// 新增类型
typedef Traits::Weighted_point_2    Weighted_point;
 
// 核心方法
// 插入加权点
Vertex_handle insert(const Weighted_point& wp, Face_handle hint = Face_handle());
 
// 点定位（基于幂距离）
Face_handle locate(const Weighted_point& wp, Face_handle hint = Face_handle());
 
// 隐藏点相关
bool is_hidden(Vertex_handle v);           // 顶点是否被隐藏
void hide_vertex(Vertex_handle v);         // 隐藏顶点
void reveal_vertex(Vertex_handle v);       // 揭示顶点
 
// 幂图相关
Point dual(Face_handle f);                 // 面的对偶点（正交中心）
Segment dual(Edge e);                      // 边的对偶线段
 
// 批量操作
template <class InputIterator>
int insert(InputIterator first, InputIterator last);

2.4 内部表示

// 顶点存储加权点
class Regular_vertex : public Delaunay_vertex {
    Weighted_point wp;
    bool hidden;
    
public:
    const Weighted_point& point() const { return wp; }
    bool is_hidden() const { return hidden; }
    void set_hidden(bool h) { hidden = h; }
};
 
// 隐藏点链表（存储在三角化数据结构中）
class Regular_triangulation_2 {
    std::list<Vertex_handle> hidden_vertices;
    // 当点被隐藏时，从三角化中移除但保留在链表中
};

3. 实现细节

3.1 加权点插入算法

// 插入加权点的主算法
Vertex_handle insert(const Weighted_point& wp) {
    // 1. 提取几何点
    Point p = wp.point();
    
    // 2. 检查是否被隐藏
    Face_handle loc = locate(wp);
    if (is_hidden_by_face(wp, loc)) {
        // 点被隐藏，加入隐藏链表
        Vertex_handle v = create_vertex(wp);
        v->set_hidden(true);
        hidden_vertices.push_back(v);
        return v;
    }
    
    // 3. 检查是否有被该点隐藏的顶点
    std::list<Vertex_handle> newly_hidden;
    find_hidden_vertices(wp, newly_hidden);
    
    // 4. 隐藏这些顶点
    for (Vertex_handle hv : newly_hidden) {
        hide_vertex(hv);
    }
    
    // 5. 插入点（使用Delaunay插入）
    Vertex_handle v = Delaunay_triangulation_2::insert(wp);
    
    // 6. 恢复正则性质（边翻转）
    restore_regular(v);
    
    return v;
}
 
// 检查点是否被面隐藏
bool is_hidden_by_face(const Weighted_point& wp, Face_handle f) {
    // 检查wp是否被f的三个顶点中的任何一个隐藏
    for (int i = 0; i < 3; ++i) {
        const Weighted_point& wp_i = f->vertex(i)->point();
        if (power_distance(wp_i, wp.point()) < 0) {
            return true;  // 被隐藏
        }
    }
    return false;
}
 
// 恢复正则性质
void restore_regular(Vertex_handle v) {
    // 类似Delaunay的边翻转，但使用power_test
    std::stack<Edge> edges_to_check;
    
    // 初始化边栈
    Face_circulator fc = incident_faces(v), done(fc);
    do {
        int i = fc->index(v);
        edges_to_check.push(Edge(fc, i));
    } while (++fc != done);
    
    // 边翻转循环
    while (!edges_to_check.empty()) {
        Edge e = edges_to_check.top();
        edges_to_check.pop();
        
        Face_handle f = e.first;
        int i = e.second;
        Face_handle n = f->neighbor(i);
        
        if (is_infinite(n)) continue;
        
        // 使用power_test代替side_of_oriented_circle
        Vertex_handle v_op = n->vertex(n->index(f));
        if (power_test(f, v_op->point()) == ON_POSITIVE_SIDE) {
            flip(f, i);
            edges_to_check.push(Edge(f, (i+1)%3));
            edges_to_check.push(Edge(n, (n->index(f)+1)%3));
        }
    }
}

3.2 幂测试 (Power Test)

// 判断点相对于正交圆的位置
Oriented_side power_test(Face_handle f, const Weighted_point& wp) {
    // 计算三角形三个加权点的正交圆
    // 然后测试wp相对于该圆的位置
    
    const Weighted_point& wp0 = f->vertex(0)->point();
    const Weighted_point& wp1 = f->vertex(1)->point();
    const Weighted_point& wp2 = f->vertex(2)->point();
    
    // 使用行列式计算
    // | x0  y0  x0²+y0²-w0  1 |
    // | x1  y1  x1²+y1²-w1  1 |
    // | x2  y2  x2²+y2²-w2  1 |
    // | xp  yp  xp²+yp²-wp  1 |
    
    return traits.power_test_2_object()(wp0, wp1, wp2, wp);
}
 
// 幂距离计算
double power_distance(const Weighted_point& wp1, const Point& p2) {
    double dx = wp1.x() - p2.x();
    double dy = wp1.y() - p2.y();
    return dx*dx + dy*dy - wp1.weight();
}

3.3 删除算法

// 删除顶点
void remove(Vertex_handle v) {
    // 1. 检查是否有被隐藏的顶点可能被揭示
    std::list<Vertex_handle> to_reveal;
    
    for (Vertex_handle hv : hidden_vertices) {
        if (would_be_revealed(hv, v)) {
            to_reveal.push_back(hv);
        }
    }
    
    // 2. 执行Delaunay删除
    Delaunay_triangulation_2::remove(v);
    
    // 3. 揭示顶点
    for (Vertex_handle hv : to_reveal) {
        reveal_vertex(hv);
        hidden_vertices.remove(hv);
        
        // 重新插入（现在可见）
        insert(hv->point());
        delete_vertex(hv);  // 删除旧的隐藏顶点
    }
}
 
// 检查删除v后hv是否会被揭示
bool would_be_revealed(Vertex_handle hv, Vertex_handle v_to_remove) {
    // 如果hv被v_to_remove隐藏，则删除后可能被揭示
    const Weighted_point& wp_hv = hv->point();
    const Weighted_point& wp_v = v_to_remove->point();
    
    return power_distance(wp_v, wp_hv.point()) < 0;
}

3.4 幂图计算

// 计算面的对偶点（正交中心）
Point power_center(Face_handle f) {
    const Weighted_point& wp0 = f->vertex(0)->point();
    const Weighted_point& wp1 = f->vertex(1)->point();
    const Weighted_point& wp2 = f->vertex(2)->point();
    
    // 计算三个加权点的正交中心
    // 这是到三个圆都有相等幂距离的点
    
    return traits.construct_power_center_2_object()(wp0, wp1, wp2);
}
 
// 计算边的对偶
Segment dual_power_edge(Edge e) {
    Face_handle f1 = e.first;
    Face_handle f2 = f1->neighbor(e.second);
    
    Point pc1 = power_center(f1);
    Point pc2 = power_center(f2);
    
    return Segment(pc1, pc2);
}

4. 完整代码示例

4.1 基础正则三角化

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Regular_triangulation_2.h>
#include <iostream>
#include <vector>
 
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Regular = CGAL::Regular_triangulation_2<K>;
 
using Weighted_point = Regular::Weighted_point;
using Point = K::Point_2;
using Vertex_handle = Regular::Vertex_handle;
using Face_handle = Regular::Face_handle;
 
int main() {
    Regular rt;
    
    // ========== 插入加权点 ==========
    std::cout << "=== 正则三角化 ===" << std::endl;
    
    // 创建加权点 (位置, 权重)
    // 权重 = 半径的平方
    std::vector<Weighted_point> weighted_points = {
        Weighted_point(Point(0, 0), 1.0),    // 圆心在(0,0)，半径1
        Weighted_point(Point(4, 0), 1.0),    // 圆心在(4,0)，半径1
        Weighted_point(Point(2, 3), 1.0),    // 圆心在(2,3)，半径1
        Weighted_point(Point(2, 1), 0.5),    // 圆心在(2,1)，半径sqrt(0.5)
        Weighted_point(Point(1, 1), 0.25),   // 圆心在(1,1)，半径0.5
        Weighted_point(Point(3, 1), 0.25),   // 圆心在(3,1)，半径0.5
    };
    
    // 批量插入
    rt.insert(weighted_points.begin(), weighted_points.end());
    
    std::cout << "插入 " << weighted_points.size() << " 个加权点" << std::endl;
    std::cout << "三角化中顶点数: " << rt.number_of_vertices() << std::endl;
    std::cout << "面数: " << rt.number_of_faces() << std::endl;
    
    // ========== 检查隐藏点 ==========
    std::cout << "\n=== 隐藏点检查 ===" << std::endl;
    
    int hidden_count = 0;
    int visible_count = 0;
    
    for (auto vit = rt.all_vertices_begin(); vit != rt.all_vertices_end(); ++vit) {
        if (rt.is_hidden(vit)) {
            ++hidden_count;
        } else {
            ++visible_count;
        }
    }
    
    std::cout << "可见点: " << visible_count << std::endl;
    std::cout << "隐藏点: " << hidden_count << std::endl;
    
    // ========== 遍历三角形 ==========
    std::cout << "\n=== 正则三角化中的三角形 ===" << std::endl;
    
    int triangle_count = 0;
    for (auto fit = rt.finite_faces_begin(); fit != rt.finite_faces_end(); ++fit) {
        std::cout << "三角形 " << ++triangle_count << ": ";
        for (int i = 0; i < 3; ++i) {
            Weighted_point wp = fit->vertex(i)->point();
            Point p = wp.point();
            std::cout << "[(" << p.x() << ", " << p.y() << "), w=" << wp.weight() << "] ";
        }
        std::cout << std::endl;
    }
    
    // ========== 幂测试验证 ==========
    std::cout << "\n=== 正则性质验证 ===" << std::endl;
    
    bool is_regular = true;
    for (auto fit = rt.finite_faces_begin(); fit != rt.finite_faces_end(); ++fit) {
        for (auto vit = rt.finite_vertices_begin(); vit != rt.finite_vertices_end(); ++vit) {
            if (fit->has_vertex(vit)) continue;
            
            auto side = rt.power_test(fit, vit->point());
            if (side == CGAL::ON_POSITIVE_SIDE) {
                std::cout << "违反正则性质!" << std::endl;
                is_regular = false;
                break;
            }
        }
    }
    std::cout << "正则验证: " << (is_regular ? "通过" : "失败") << std::endl;
    
    // ========== 点定位 ==========
    std::cout << "\n=== 点定位 ===" << std::endl;
    
    Weighted_point query(Point(2, 1.5), 0.1);
    Face_handle loc = rt.locate(query);
    
    if (!rt.is_infinite(loc)) {
        std::cout << "查询点位于三角形内，顶点为:" << std::endl;
        for (int i = 0; i < 3; ++i) {
            Weighted_point wp = loc->vertex(i)->point();
            Point p = wp.point();
            std::cout << "  [(" << p.x() << ", " << p.y() << "), w=" << wp.weight() << "]" << std::endl;
        }
    }
    
    return 0;
}

4.2 幂图计算

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Regular_triangulation_2.h>
#include <iostream>
#include <vector>
#include <fstream>
 
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Regular = CGAL::Regular_triangulation_2<K>;
 
using Weighted_point = Regular::Weighted_point;
using Point = K::Point_2;
using Segment = K::Segment_2;
using Face_handle = Regular::Face_handle;
using Edge = Regular::Edge;
 
int main() {
    Regular rt;
    
    // 创建加权点（模拟不同大小的圆）
    std::vector<Weighted_point> points;
    
    // 中心大圆
    points.push_back(Weighted_point(Point(5, 5), 4.0));
    
    // 周围小圆
    for (int i = 0; i < 8; ++i) {
        double angle = 2 * M_PI * i / 8;
        double x = 5 + 3 * cos(angle);
        double y = 5 + 3 * sin(angle);
        points.push_back(Weighted_point(Point(x, y), 0.5));
    }
    
    rt.insert(points.begin(), points.end());
    
    std::cout << "=== 幂图计算 ===" << std::endl;
    std::cout << "加权点数: " << points.size() << std::endl;
    std::cout << "三角化顶点数: " << rt.number_of_vertices() << std::endl;
    
    // ========== 计算幂图顶点（正交中心） ==========
    std::cout << "\n=== 幂图顶点 ===" << std::endl;
    
    std::vector<Point> power_vertices;
    for (auto fit = rt.finite_faces_begin(); fit != rt.finite_faces_end(); ++fit) {
        Point pc = rt.dual(fit);  // 正交中心
        power_vertices.push_back(pc);
        std::cout << "正交中心: (" << pc.x() << ", " << pc.y() << ")" << std::endl;
    }
    
    // ========== 计算幂图边 ==========
    std::cout << "\n=== 幂图边 ===" << std::endl;
    
    std::vector<Segment> power_edges;
    for (auto eit = rt.finite_edges_begin(); eit != rt.finite_edges_end(); ++eit) {
        Face_handle f1 = eit->first;
        Face_handle f2 = f1->neighbor(eit->second);
        
        if (!rt.is_infinite(f2)) {
            Segment seg = rt.dual(*eit);
            power_edges.push_back(seg);
            std::cout << "幂边: (" << seg.source().x() << ", " << seg.source().y() 
                      << ") -> (" << seg.target().x() << ", " << seg.target().y() << ")" << std::endl;
        }
    }
    
    // ========== 输出到SVG ==========
    std::ofstream svg("power_diagram.svg");
    svg << "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n";
    svg << "<svg xmlns=\"http://www.w3.org/2000/svg\" width=\"600\" height=\"600\">\n";
    
    // 绘制加权圆（输入圆）
    svg << "<g fill=\"none\" stroke=\"red\" stroke-width=\"1\">\n";
    for (const auto& wp : points) {
        Point p = wp.point();
        double r = std::sqrt(wp.weight());
        svg << "<circle cx=\"" << p.x() * 50 + 50 
            << "\" cy=\"" << 550 - p.y() * 50
            << "\" r=\"" << r * 50 << "\"/>\n";
    }
    svg << "</g>\n";
    
    // 绘制正则三角化（灰色虚线）
    svg << "<g stroke=\"gray\" stroke-width=\"1\" stroke-dasharray=\"5,5\">\n";
    for (auto eit = rt.finite_edges_begin(); eit != rt.finite_edges_end(); ++eit) {
        Segment seg = rt.segment(*eit);
        svg << "<line x1=\"" << seg.source().x() * 50 + 50
            << "\" y1=\"" << 550 - seg.source().y() * 50
            << "\" x2=\"" << seg.target().x() * 50 + 50
            << "\" y2=\"" << 550 - seg.target().y() * 50 << "\"/>\n";
    }
    svg << "</g>\n";
    
    // 绘制幂图（蓝色）
    svg << "<g stroke=\"blue\" stroke-width=\"2\">\n";
    for (const auto& seg : power_edges) {
        svg << "<line x1=\"" << seg.source().x() * 50 + 50
            << "\" y1=\"" << 550 - seg.source().y() * 50
            << "\" x2=\"" << seg.target().x() * 50 + 50
            << "\" y2=\"" << 550 - seg.target().y() * 50 << "\"/>\n";
    }
    svg << "</g>\n";
    
    // 绘制圆心（红色点）
    svg << "<g fill=\"red\">\n";
    for (const auto& wp : points) {
        Point p = wp.point();
        svg << "<circle cx=\"" << p.x() * 50 + 50
            << "\" cy=\"" << 550 - p.y() * 50 << "\" r=\"3\"/>\n";
    }
    svg << "</g>\n";
    
    svg << "</svg>\n";
    svg.close();
    
    std::cout << "\n幂图已保存到 power_diagram.svg" << std::endl;
    
    return 0;
}

4.3 隐藏点演示

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Regular_triangulation_2.h>
#include <iostream>
#include <vector>
 
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Regular = CGAL::Regular_triangulation_2<K>;
 
using Weighted_point = Regular::Weighted_point;
using Point = K::Point_2;
using Vertex_handle = Regular::Vertex_handle;
 
int main() {
    Regular rt;
    
    std::cout << "=== 隐藏点演示 ===" << std::endl;
    
    // 首先插入一个大圆
    Weighted_point big_circle(Point(5, 5), 9.0);  // 半径3
    Vertex_handle v_big = rt.insert(big_circle);
    
    std::cout << "插入大圆 [(5,5), r=3]" << std::endl;
    std::cout << "顶点数: " << rt.number_of_vertices() << std::endl;
    
    // 插入一个小圆（在大圆内部）
    Weighted_point small_circle(Point(5, 5), 0.25);  // 半径0.5，同圆心
    Vertex_handle v_small = rt.insert(small_circle);
    
    std::cout << "\n插入小圆 [(5,5), r=0.5]" << std::endl;
    std::cout << "小圆是否隐藏: " << (rt.is_hidden(v_small) ? "是" : "否") << std::endl;
    std::cout << "顶点数: " << rt.number_of_vertices() << std::endl;
    
    // 插入另一个圆（部分与大圆重叠）
    Weighted_point overlap_circle(Point(6, 5), 1.0);  // 半径1
    Vertex_handle v_overlap = rt.insert(overlap_circle);
    
    std::cout << "\n插入重叠圆 [(6,5), r=1]" << std::endl;
    std::cout << "重叠圆是否隐藏: " << (rt.is_hidden(v_overlap) ? "是" : "否") << std::endl;
    std::cout << "顶点数: " << rt.number_of_vertices() << std::endl;
    
    // 显示所有可见点
    std::cout << "\n=== 可见点 ===" << std::endl;
    for (auto vit = rt.finite_vertices_begin(); vit != rt.finite_vertices_end(); ++vit) {
        Weighted_point wp = vit->point();
        Point p = wp.point();
        std::cout << "[(" << p.x() << ", " << p.y() << "), r=" 
                  << std::sqrt(wp.weight()) << "]" << std::endl;
    }
    
    // 演示删除后揭示隐藏点
    std::cout << "\n=== 删除大圆 ===" << std::endl;
    rt.remove(v_big);
    
    std::cout << "删除后顶点数: " << rt.number_of_vertices() << std::endl;
    
    // 检查之前隐藏的点是否被揭示
    std::cout << "\n删除后的可见点:" << std::endl;
    for (auto vit = rt.finite_vertices_begin(); vit != rt.finite_vertices_end(); ++vit) {
        Weighted_point wp = vit->point();
        Point p = wp.point();
        std::cout << "[(" << p.x() << ", " << p.y() << "), r=" 
                  << std::sqrt(wp.weight()) << "]" << std::endl;
    }
    
    return 0;
}

5. 复杂度分析

5.1 时间复杂度

操作平均复杂度最坏复杂度说明
insert(wp)O(log n)O(n)加权点插入
insert(first, last)O(n log n)O(n²)批量插入
remove(v)O(log n)O(n)顶点删除（可能揭示隐藏点）
locate(wp)O(log n)O(n)点定位
is_hidden(v)O(1)O(1)常数时间
power_testO(1)O(1)正交圆测试

5.2 空间复杂度

存储项目空间
基础三角化O(n)
隐藏点链表O(h)，h为隐藏点数
权重存储O(n)
总计O(n + h)

5.3 与Delaunay的比较

特性DelaunayRegular
输入普通点加权点
插入O(log n)O(log n)
隐藏点无有
删除复杂度O(log n)O(log n + r)，r为揭示点数
空间开销O(n)O(n + h)
应用场景一般三角化带权重的邻近查询、球体填充

6. 应用场景

6.1 分子生物学

// 计算分子表面的溶剂可及表面积
// 原子作为加权点（位置+范德华半径）
// 正则三角化计算分子表面

6.2 球体填充与堆积

// 不同大小球体的空间填充
// 球心作为加权点，半径平方作为权重
// 正则三角化给出邻近关系

6.3 各向异性网格生成

// 根据局部特征尺寸调整权重
// 生成自适应网格

6.4 经济学与博弈论

// 权力指数计算
// 市场区域划分

7. 总结

Regular_triangulation_2是Delaunay三角化在加权点上的自然推广：

  1. 统一的框架：当所有权重为0时退化为Delaunay三角化
  2. 隐藏点机制：自动处理被完全包含的加权点
  3. 幂图对偶：支持加权点的邻近查询和空间分解
  4. 广泛的应用：从分子生物学到经济学

正则三角化展示了CGAL三角化框架的灵活性，通过模板参数和继承机制，在保持接口一致性的同时支持复杂的几何概念。

Graph View

Backlinks