基于Python的CGAL裁剪Voronoi图至边界框及相关问题咨询

背景 #

我正在用Python结合CGAL绘制Voronoi图，需要计算Voronoi图被限制在一个边界框内时，特定站点周围的边集合。目前已经实现了获取特定站点周围有界边的逆时针遍历，但无法处理无界边的情况（无界边的源或目标为无穷远，用None避免内核崩溃）。

以下是我编写的代码：

import numpy as np
from CGAL import CGAL_Kernel
from CGAL.CGAL_Voronoi_diagram_2 import Voronoi_diagram_2
from CGAL.CGAL_Triangulation_2 import Delaunay_triangulation_2

sites = np.array([[-3, 2], [-5, -3], [2,5], [0,5], [0,0]], dtype=float)
print("sites = ", sites,end="\n\n")
sites_cgal = [CGAL_Kernel.Point_2(x, y) for x, y in sites]
dt = Delaunay_triangulation_2()
dt.insert(sites_cgal)
vd = Voronoi_diagram_2(dt)

for f in vd.faces():
    e = f.halfedge()
    d = f.dual()
    print('site = ', d.point())
    s = [e.source().point().x(),e.source().point().y()] if e.has_source() else [None,None]
    t = [e.target().point().x(),e.target().point().y()] if e.has_target() else [None,None]
    print("edge : source = ",s," target = ", t)
    while True:
        e = e.next()
        sn = [e.source().point().x(),e.source().point().y()] if e.has_source() else [None,None]
        tn = [e.target().point().x(),e.target().point().y()] if e.has_target() else [None,None]
        if [s,t] == [sn,tn]:
            break
        else:
            print("edge : source = ", sn," target = ", tn)
    
    print("\n")

需要解决的问题 #

1. 裁剪源/目标为无穷远的无界边，找出每个站点对应的无界边与边界框相交后的边；
2. 了解除源、目标外，是否有其他元素可将无界边转换为射线、线段等形式；
3. 咨询是否可使用from CGAL.CGAL_Voronoi_diagram_2 import crop_voronoi_edge函数裁剪Voronoi图，若可以，具体如何使用？

解决方案 #

问题1：裁剪无界边与边界框的相交边 #

无界Voronoi边本质是平分线射线，对应Delaunay三角剖分中与无穷远点相连的边。裁剪这类边可按以下步骤操作：

• 定义边界框（用CGAL_Kernel.Rectangle_2或轴对齐矩形的四条边表示）；
• 获取无界边的支撑线：通过边对应的Delaunay边拿到两个相邻站点，计算平分线的中点和方向，生成射线；
• 求射线与边界框的交点，将无界边转换为边界框内的线段。

示例代码片段（边界框设为x∈[-6,3], y∈[-4,6]）：

from CGAL.CGAL_Kernel import Line_2, Ray_2, Point_2, Segment_2

# 定义边界框的四条边
bbox_edges = [
    Segment_2(Point_2(-6, -4), Point_2(3, -4)),  # 底边
    Segment_2(Point_2(3, -4), Point_2(3, 6)),    # 右边
    Segment_2(Point_2(3, 6), Point_2(-6, 6)),    # 顶边
    Segment_2(Point_2(-6, 6), Point_2(-6, -4))   # 左边
]

def crop_unbounded_edge(v_edge, bbox_edges):
    # 获取Voronoi边对应的两个站点
    delaunay_edge = v_edge.dual()
    p1 = delaunay_edge.first().point()
    p2 = delaunay_edge.second().point()
    # 计算平分线的中垂线（支撑线）
    midpoint = Point_2((p1.x()+p2.x())/2, (p1.y()+p2.y())/2)
    direction = Point_2(p1.y()-p2.y(), p2.x()-p1.x())  # 垂直于p1-p2的方向
    ray = Ray_2(midpoint, direction)
    
    # 找射线与边界框的第一个交点
    intersection_point = None
    for seg in bbox_edges:
        result = ray.intersection(seg)
        if result is not None:
            if isinstance(result, Point_2):
                intersection_point = result
                break
            # 若交点是线段（射线与边界边重合，取远端点）
            elif isinstance(result, Segment_2):
                dist1 = ray.source().squared_distance(result.source())
                dist2 = ray.source().squared_distance(result.target())
                intersection_point = result.target() if dist2 > dist1 else result.source()
                break
    # 返回裁剪后的线段
    if intersection_point:
        return Segment_2(midpoint, intersection_point)
    return None

问题2：无界边转射线/线段的其他方式 #

CGAL的Voronoi边自带接口可直接转换为射线：

• 用v_edge.is_unbounded()判断边是否无界；
• 用v_edge.line()获取边的支撑线；
• 结合Voronoi面的逆时针方向，确定射线的起点（存在的Voronoi顶点）和延伸方向：
  • 若边的源为无穷远，射线起点为目标顶点，方向指向无穷；
  • 若边的目标为无穷远，射线起点为源顶点，方向指向无穷。

示例代码：

for f in vd.faces():
    e = f.halfedge()
    start_e = e
    while True:
        if e.is_unbounded():
            # 确定射线起点
            start = e.target().point() if not e.has_source() else e.source().point()
            # 获取支撑线并确定方向
            line = e.line()
            direction = line.direction() if f.is_ccw_around_point(start) else line.opposite().direction()
            ray = Ray_2(start, direction)
            print(f"Unbounded ray: {start} -> {direction}")
        e = e.next()
        if e == start_e:
            break

问题3：使用crop_voronoi_edge函数 #

crop_voronoi_edge可以直接裁剪Voronoi边（支持有界/无界边），返回边界框内的线段（无交点则返回None）。使用步骤如下：

1. 导入函数：from CGAL.CGAL_Voronoi_diagram_2 import crop_voronoi_edge
2. 定义边界框为CGAL_Kernel.Rectangle_2对象；
3. 遍历Voronoi边并调用函数裁剪。

示例代码：

from CGAL.CGAL_Kernel import Rectangle_2, Point_2

# 定义边界框
bbox = Rectangle_2(Point_2(-6, -4), Point_2(3, 6))

for f in vd.faces():
    e = f.halfedge()
    site = f.dual().point()
    print(f"Site: {site}")
    start_e = e
    while True:
        if e.is_bounded():
            # 有界边直接输出
            seg = Segment_2(e.source().point(), e.target().point())
            print(f"Bounded edge: {seg.source()} -> {seg.target()}")
        else:
            # 无界边裁剪
            cropped_seg = crop_voronoi_edge(e, bbox)
            if cropped_seg:
                print(f"Cropped unbounded edge: {cropped_seg.source()} -> {cropped_seg.target()}")
        e = e.next()
        if e == start_e:
            break
    print("\n")

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