如何为含贝塞尔曲线的matplotlib.PathPatch计算给定原点的最小覆盖圆？

数据老陈

2026-7-24

求解覆盖Matplotlib PathPatch的最小半径圆（含贝塞尔曲线场景）

问题描述

针对任意Matplotlib PathPatch实例与给定原点，需要找到能完全覆盖该图形的最小半径圆：

对于仅由直线构成的路径，解法很直接：计算所有路径顶点到原点的距离，取最大值即为所需半径。
但对于包含贝塞尔曲线的路径，路径顶点是控制点，可能位于图形外部，直接用控制点计算出的圆半径会过大。需要找到曲线段上距离原点最远的点，才能得到真正的最小覆盖半径。

效果示意图

核心思路

覆盖圆的最小半径，等于路径上所有点到原点的最大距离。对于贝塞尔曲线段，这个最大距离可能出现在曲线内部而非控制点上，因此需要对每个贝塞尔线段进行极值求解，再结合路径顶点的距离，取最大值。

具体实现步骤

1. 解析路径线段

Matplotlib的Path对象通过vertices和codes定义路径，需根据代码拆分出不同类型的线段：

Path.MOVETO(1)：仅移动起点，不构成线段
Path.LINETO(2)：直线段，连接当前点到下一个顶点
Path.CURVE3(3)：二次贝塞尔曲线，需当前点+2个控制点
Path.CURVE4(4)：三次贝塞尔曲线，需当前点+3个控制点
Path.CLOSEPOLY(79)：闭合路径，连接当前点到起点

2. 计算各类线段的最大距离

为简化计算，我们求距离平方的极值（与距离的极值点一致，避免开根号运算）：

直线段：最大距离必然出现在端点，直接计算两个端点的距离取最大值。
二次贝塞尔曲线：对距离平方的参数方程求导，解出[0,1]范围内的极值点参数t，计算这些点和端点的距离，取最大值。
三次贝塞尔曲线：对距离平方的参数方程求导得到三次方程，解出[0,1]范围内的实数根，计算对应点的距离，再和端点距离比较取最大值。

3. 整合所有线段的最大距离

遍历路径所有线段，收集每个线段的最大距离，最终的最大值就是最小覆盖圆的半径。

完整代码实现

#!/usr/bin/env python
"""
Find the circle with minimum radius that given an xy origin fully covers a matplotlib.patch.PathPatch instance.
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch

def distance_sq(point, origin):
    """计算点到原点的距离平方"""
    return np.sum((point - origin)**2)

def max_distance_line(p0, p1, origin):
    """计算直线段上到原点的最大距离"""
    d0 = np.sqrt(distance_sq(p0, origin))
    d1 = np.sqrt(distance_sq(p1, origin))
    return max(d0, d1)

def max_distance_quadratic_bezier(p0, p1, p2, origin):
    """计算二次贝塞尔曲线上到原点的最大距离"""
    p0_rel = p0 - origin
    p1_rel = p1 - origin
    p2_rel = p2 - origin

    # 推导距离平方导数为0的二次方程系数
    coeff_A = np.dot(p2_rel - p0_rel, p2_rel - 2*p1_rel + p0_rel)
    coeff_B = 2 * np.dot(p1_rel - p0_rel, p2_rel - 2*p1_rel + p0_rel)
    coeff_C = np.dot(p1_rel - p0_rel, p1_rel - p0_rel) + np.dot(p0_rel, p2_rel - p1_rel)

    candidates_t = []
    discriminant = coeff_B**2 - 4*coeff_A*coeff_C
    if discriminant >= 0:
        sqrt_disc = np.sqrt(discriminant)
        if coeff_A != 0:
            t1 = (-coeff_B - sqrt_disc) / (2*coeff_A)
            t2 = (-coeff_B + sqrt_disc) / (2*coeff_A)
            for t in [t1, t2]:
                if 0 <= t <= 1:
                    candidates_t.append(t)
        else:
            if coeff_B != 0:
                t = -coeff_C / coeff_B
                if 0 <= t <= 1:
                    candidates_t.append(t)
    
    candidates_t.extend([0, 1])

    max_dist = 0
    for t in candidates_t:
        pt = (1-t)**2 * p0 + 2*t*(1-t)*p1 + t**2 * p2
        dist = np.sqrt(distance_sq(pt, origin))
        if dist > max_dist:
            max_dist = dist
    return max_dist

def max_distance_cubic_bezier(p0, p1, p2, p3, origin):
    """计算三次贝塞尔曲线上到原点的最大距离"""
    p0_rel = p0 - origin
    p1_rel = p1 - origin
    p2_rel = p2 - origin
    p3_rel = p3 - origin

    # 推导距离平方导数为0的三次方程系数
    a, b, c, d = p0_rel, p1_rel, p2_rel, p3_rel
    coeffs = np.poly1d([
        np.dot(d - a, d - 3*c + 3*b - a),
        3*np.dot(d - a, c - 2*b + a) + 3*np.dot(c - a, d - 3*c + 3*b - a),
        3*np.dot(d - a, b - a) + 6*np.dot(c - a, c - 2*b + a) + 3*np.dot(b - a, d - 3*c + 3*b - a),
        3*np.dot(c - a, b - a) + 3*np.dot(b - a, c - 2*b + a)
    ]).coeffs

    roots = np.roots(coeffs)
    candidates_t = []
    for root in roots:
        if np.isreal(root):
            t = np.real(root)
            if 0 <= t <= 1:
                candidates_t.append(t)
    
    candidates_t.extend([0, 1])

    max_dist = 0
    for t in candidates_t:
        pt = (1-t)**3 * p0 + 3*t*(1-t)**2 * p1 + 3*t**2*(1-t)*p2 + t**3 * p3
        dist = np.sqrt(distance_sq(pt, origin))
        if dist > max_dist:
            max_dist = dist
    return max_dist

def find_min_covering_radius(path_patch, origin):
    """计算覆盖PathPatch的最小圆半径"""
    path = path_patch.get_path()
    vertices = path.vertices
    codes = path.codes
    origin = np.asarray(origin)

    max_radius = 0
    n = len(vertices)
    i = 0
    while i < n:
        code = codes[i]
        if code == Path.MOVETO:
            dist = np.sqrt(distance_sq(vertices[i], origin))
            max_radius = max(max_radius, dist)
            i += 1
        elif code == Path.LINETO:
            p0 = vertices[i-1]
            p1 = vertices[i]
            dist = max_distance_line(p0, p1, origin)
            max_radius = max(max_radius, dist)
            i += 1
        elif code == Path.CURVE3:
            p0 = vertices[i-1]
            p1 = vertices[i]
            p2 = vertices[i+1]
            dist = max_distance_quadratic_bezier(p0, p1, p2, origin)
            max_radius = max(max_radius, dist)
            i += 2
        elif code == Path.CURVE4:
            p0 = vertices[i-1]
            p1 = vertices[i]
            p2 = vertices[i+1]
            p3 = vertices[i+2]
            dist = max_distance_cubic_bezier(p0, p1, p2, p3, origin)
            max_radius = max(max_radius, dist)
            i += 3
        elif code == Path.CLOSEPOLY:
            p0 = vertices[i-1]
            p1 = vertices[0]
            dist = max_distance_line(p0, p1, origin)
            max_radius = max(max_radius, dist)
            i += 1
        else:
            i += 1
    return max_radius

# ---------------------- 测试示例 ----------------------
fig, ax = plt.subplots()

# 构造PathPatch
vertices = np.array([[ 0.44833333, -2.75444444],
                     [-0.78166667, -1.28444444],
                     [-2.88166667,  1.81555556],
                     [-0.75666667,  1.81555556],
                     [-0.28166667,  0.96555556],
                     [ 1.06833333,  3.01555556],
                     [ 1.86833333, -0.13444444],
                     [ 0.86833333, -0.68444444],
                     [ 0.44833333, -2.75444444]])
codes = (1, 4, 4, 4, 2, 4, 4, 4, 79)
path_patch = PathPatch(Path(vertices, codes), color='red')
ax.add_artist(path_patch)

# 绘制控制点
ax.scatter(*vertices.T, c='black', marker='x')

# 原点
origin = np.array((0, 0))

# 原方法计算的半径（仅用控制点）
deltas = vertices - origin[np.newaxis, :]
distances = np.linalg.norm(deltas, axis=-1)
old_radius = np.max(distances)
ax.add_artist(plt.Circle(origin, old_radius, alpha=0.1, color='blue', label='原方法半径'))

# 新方法计算的最小半径
min_radius = find_min_covering_radius(path_patch, origin)
ax.add_artist(plt.Circle(origin, min_radius, alpha=0.2, color='green', label='最小覆盖半径'))

# 设置坐标轴
ax.axis([-4, 4, -4, 4])
ax.set_aspect("equal")
plt.legend()
plt.show()

print(f"原方法半径：{old_radius:.3f}")
print(f"最小覆盖半径：{min_radius:.3f}")

代码说明

距离平方计算：用距离平方代替距离，避免开根号操作，提升计算效率且极值点一致。
贝塞尔曲线极值求解：通过对距离平方的参数方程求导，找到极值点对应的参数t，再计算该点的距离。对于三次贝塞尔曲线，使用numpy求解三次方程的根，筛选出有效范围内的实数根。
路径遍历：根据Path的codes逐个解析线段类型，分别处理直线、二次/三次贝塞尔曲线，收集所有线段的最大距离。

效果对比

运行代码后，蓝色圆是原方法（仅用控制点）得到的过大半径，绿色圆是真正的最小覆盖圆，能恰好包围整个PathPatch图形。

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

如何为含贝塞尔曲线的matplotlib.PathPatch计算给定原点的最小覆盖圆？

问题描述 #

核心思路 #

具体实现步骤 #

1. 解析路径线段 #

2. 计算各类线段的最大距离 #

3. 整合所有线段的最大距离 #

完整代码实现 #

代码说明 #

效果对比 #