基于Matplotlib开发Python子空间绘制函数的技术咨询

先搞定线性无关性检查 #

把输入的向量转成矩阵，用numpy.linalg.matrix_rank()计算秩，只要秩等于向量个数，就说明线性无关，可继续执行；否则直接抛错提示输入无效。注意所有向量必须是3维的，适配3D绘图场景。

轴边界裁剪的核心逻辑 #

先定义轴的边界范围，比如x_min, x_max, y_min, y_max, z_min, z_max，所有绘制的线、面、体都必须限制在这个范围内。

1. 单向量：绘制裁剪后的直线 #

• 直线用参数方程表示：r(t) = t * v（默认向量起点在原点，若有偏移可自行调整）
• 求解t的有效范围：将x、y、z分别代入边界，算出对应的t值，取这些t值的交集，得到t_start和t_end
• 计算两个端点t_start*v和t_end*v，直接用ax.plot()绘制线段即可，无需复杂集合操作。

2. 两个向量：绘制裁剪后的平面 #

这是核心难点，需要找到平面与轴边界6个面的交线，再将交线的有效端点整理成有序多边形：

• 先推导平面方程：两个向量v1、v2的叉乘即为平面法向量n，平面方程为n[0]x + n[1]y + n[2]z = 0（原点在平面上）
• 遍历6个轴边界平面（如x=x_min、x=x_max等），求解当前平面与每个边界平面的交线：
  • 若两个平面平行（法向量点积接近0），直接跳过；否则计算交线的方向向量（两个法向量的叉乘）和一个点（联立方程求解）
• 对每条交线，复用单向量的裁剪逻辑，得到边界内的线段端点
• 收集所有有效端点后，自动排序顶点：
  • 用scipy.spatial.ConvexHull求凸包是最可靠的方法，凸包输出的顶点顺序天然是能构成闭合多边形的正确顺序，无需手动计算极角，避免出错
• 最后将排序后的顶点列表传给Poly3DCollection，记得把第一个点添加到列表末尾，让多边形闭合。

3. 三个向量：绘制轴边界的立方体 #

这部分本质是绘制3D边界框，直接用Poly3DCollection绘制6个面即可：

• 先定义8个顶点：比如[(x_min,y_min,z_min), (x_min,y_min,z_max), ..., (x_max,y_max,z_max)]
• 每个面对应4个顶点，注意所有面的顶点顺序要统一（比如从外部看均为顺时针），否则会出现面显示异常
• 将6个面的顶点列表传入Poly3DCollection，设置颜色和透明度即可。

关键细节提醒 #

• 交点计算通用方法：不管是直线还是平面与边界的交点，都可通过联立方程求解，比如将直线参数方程代入边界平面，解出参数t后，验证对应点是否在其他轴的边界内
• 顶点排序避坑：不要手动排序顶点，用凸包算法自动处理，尤其是平面与多个边界相交时，凸包能自动筛选出外围顶点，保证多边形为凸且顺序正确
• Poly3DCollection注意事项：每个面的顶点必须闭合，且所有面的顶点顺序要一致（全顺时针或全逆时针），否则会出现面的正反面显示混乱；可通过edgecolors设置边框颜色，让结构更清晰。

代码框架示例 #

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from scipy.spatial import ConvexHull

def check_linear_independence(vectors):
    mat = np.array(vectors).T
    rank = np.linalg.matrix_rank(mat)
    return rank == len(vectors)

def clip_line(v, bounds):
    x_min, x_max, y_min, y_max, z_min, z_max = bounds
    t_values = []
    # 计算所有可能的t值
    if abs(v[0]) > 1e-9:
        t_values.extend([x_min / v[0], x_max / v[0]])
    if abs(v[1]) > 1e-9:
        t_values.extend([y_min / v[1], y_max / v[1]])
    if abs(v[2]) > 1e-9:
        t_values.extend([z_min / v[2], z_max / v[2]])
    t_values = np.array(t_values)
    # 筛选出在所有轴边界内的t值（允许微小误差）
    t_valid = t_values[
        (t_values * v[0] >= x_min - 1e-6) & (t_values * v[0] <= x_max + 1e-6) &
        (t_values * v[1] >= y_min - 1e-6) & (t_values * v[1] <= y_max + 1e-6) &
        (t_values * v[2] >= z_min - 1e-6) & (t_values * v[2] <= z_max + 1e-6)
    ]
    if len(t_valid) < 2:
        return None
    t_start, t_end = t_valid.min(), t_valid.max()
    return np.array([t_start * v, t_end * v])

def clip_plane(v1, v2, bounds):
    x_min, x_max, y_min, y_max, z_min, z_max = bounds
    n = np.cross(v1, v2)
    intersections = []
    # 遍历6个边界平面
    for axis, val in [
        ('x', x_min), ('x', x_max),
        ('y', y_min), ('y', y_max),
        ('z', z_min), ('z', z_max)
    ]:
        # 求当前边界平面与目标平面的交线方向向量
        plane_normal = np.array([1,0,0]) if axis == 'x' else np.array([0,1,0]) if axis == 'y' else np.array([0,0,1])
        dir_vec = np.cross(n, plane_normal)
        if np.linalg.norm(dir_vec) < 1e-9:
            continue  # 平面平行，无交线
        # 找交线上的一个点
        if axis == 'x':
            y = 0
            z = (-n[0] * val) / n[2] if abs(n[2]) > 1e-9 else 0
            if abs(n[1]) > 1e-9 and abs(n[2]) < 1e-9:
                z = 0
                y = (-n[0] * val) / n[1]
            point = np.array([val, y, z])
        elif axis == 'y':
            x = 0
            z = (-n[1] * val) / n[2] if abs(n[2]) > 1e-9 else 0
            if abs(n[0]) > 1e-9 and abs(n[2]) < 1e-9:
                z = 0
                x = (-n[1] * val) / n[0]
            point = np.array([x, val, z])
        else:
            x = 0
            y = (-n[2] * val) / n[1] if abs(n[1]) > 1e-9 else 0
            if abs(n[0]) > 1e-9 and abs(n[1]) < 1e-9:
                y = 0
                x = (-n[2] * val) / n[0]
            point = np.array([x, y, val])
        # 裁剪交线，得到有效端点
        line_seg = clip_line(dir_vec, bounds)
        if line_seg is not None:
            # 验证端点是否在目标平面上（允许微小误差）
            for p in line_seg:
                if abs(np.dot(n, p)) < 1e-6:
                    intersections.append(p)
    if len(intersections) < 3:
        return None
    # 用凸包获取有序的多边形顶点
    hull = ConvexHull(intersections)
    ordered_verts = intersections[hull.vertices]
    # 闭合多边形
    ordered_verts = np.vstack([ordered_verts, ordered_verts[0]])
    return ordered_verts

def get_bounds_verts(bounds):
    x_min, x_max, y_min, y_max, z_min, z_max = bounds
    # 定义8个顶点
    v0 = (x_min, y_min, z_min)
    v1 = (x_min, y_min, z_max)
    v2 = (x_min, y_max, z_min)
    v3 = (x_min, y_max, z_max)
    v4 = (x_max, y_min, z_min)
    v5 = (x_max, y_min, z_max)
    v6 = (x_max, y_max, z_min)
    v7 = (x_max, y_max, z_max)
    # 6个面的顶点列表（顺序统一为顺时针，从外部看）
    return [
        [v0, v1, v3, v2],  # x=x_min面
        [v4, v5, v7, v6],  # x=x_max面
        [v0, v4, v6, v2],  # z=z_min面
        [v1, v5, v7, v3],  # z=z_max面
        [v0, v4, v5, v1],  # y=y_min面
        [v2, v6, v7, v3],  # y=y_max面
    ]

def plot_subspace(vectors, bounds=(-1,1,-1,1,-1,1)):
    if not check_linear_independence(vectors):
        raise ValueError("输入向量必须线性无关")
    fig = plt.figure(figsize=(8,6))
    ax = fig.add_subplot(111, projection='3d')
    ax.set_xlim(bounds[0], bounds[1])
    ax.set_ylim(bounds[2], bounds[3])
    ax.set_zlim(bounds[4], bounds[5])
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    
    if len(vectors) == 1:
        line_points = clip_line(vectors[0], bounds)
        if line_points is not None:
            ax.plot(line_points[:,0], line_points[:,1], line_points[:,2], color='#1f77b4', linewidth=2)
    elif len(vectors) == 2:
        plane_verts = clip_plane(vectors[0], vectors[1], bounds)
        if plane_verts is not None:
            poly = Poly3DCollection([plane_verts], facecolors='#ff7f0e', edgecolors='#1f77b4', alpha=0.5)
            ax.add_collection3d(poly)
    elif len(vectors) == 3:
        box_verts = get_bounds_verts(bounds)
        poly = Poly3DCollection(box_verts, facecolors=['#1f77b4','#1f77b4','#ff7f0e','#ff7f0e','#2ca02c','#2ca02c'], 
                                edgecolors='black', alpha=0.3)
        ax.add_collection3d(poly)
    plt.tight_layout()
    plt.show()

# 测试用例
plot_subspace([np.array([1,2,3])], bounds=(-2,2,-2,2,-2,2))
plot_subspace([np.array([1,0,0]), np.array([0,1,0])], bounds=(-1,1,-1,1,-1,1))
plot_subspace([np.array([1,0,0]), np.array([0,1,0]), np.array([0,0,1])], bounds=(-1,1,-1,1,-1,1))

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