如何在3D轨迹图上绘制随点移动且垂直于蓝色坐标轴的圆

问题背景 #

现有代码可生成沿轨迹移动的点，以及随点旋转的三个正交坐标轴（红色切向量T、蓝色副法向量B、绿色法向量N）。需要在原3D图中添加垂直于蓝色B向量的圆，且该圆能随轨迹点移动（随theta角和h变量）自动调整倾斜角度。

实现思路 #

1. 圆所在平面的法向量就是当前点的B向量，圆心为轨迹上的当前点(x[t-1], y[t-1], z[t-1])
2. 利用T和N向量作为圆所在平面的两个正交单位基（T、N、B构成正交标架），通过参数方程生成圆的所有点
3. 将生成的圆的基础点平移到当前圆心位置，即可得到随点倾斜的圆

修改后的完整代码 #

import numpy as np
import matplotlib.pyplot as plt

r = 18 
h = 1.7
num_of_steps = 1000
emp = 3
time = np.arange(0, 100, 1)
phi = []
theta = []
Amp = np.pi/8
fphi = 4
ftheta = 2
pics = []

for j in time:
    kampas = np.radians(2*np.pi*fphi*j)
    kitaskampas = Amp*(np.sin(np.radians(2*np.pi*ftheta*j)))
    phi.append(kampas)
    theta.append(kitaskampas)
# 生成随时间变化的phi和theta角度列表

x = r * np.cos(phi)
y = r * np.sin(phi) * np.cos(theta) - h * np.sin(theta)
z = r * np.sin(phi) * np.sin(theta) + h * np.cos(theta)
# 计算轨迹点的x、y、z坐标
points = np.transpose([x, y, z])

n = len(points)

# 计算轨迹的一阶和二阶导数
dX = np.apply_along_axis(np.gradient, axis=0, arr=points)
ddX = np.apply_along_axis(np.gradient, axis=0, arr=dX)

# 归一化切向量T（增加零向量判断避免除零错误）
f = lambda m: m / np.linalg.norm(m) if np.linalg.norm(m) !=0 else m
T = np.apply_along_axis(f, axis=1, arr=dX)

# 计算并归一化副法向量B
B = np.cross(dX, ddX)
B = np.apply_along_axis(f, axis=1, arr=B)

# 计算法向量N
N = np.cross(B, T)

# 提前生成圆的基础参数点，避免循环内重复计算
circle_radius = 5  # 可自行调整圆的半径
theta_circle = np.linspace(0, 2*np.pi, 50)
circle_base = np.array([circle_radius*np.cos(theta_circle), circle_radius*np.sin(theta_circle), np.zeros_like(theta_circle)])

for t in range(1, len(theta)):  # 从1开始避免t-1索引越界
    fig = plt.figure('Parametrinai blynai')
    ax = fig.add_subplot(111, projection='3d')
    ax.clear()
    ax.plot(x[1:t], y[1:t], z[1:t], '-r', linewidth=3)
    ax.azim = 0
    ax.dist = 10
    ax.elev = 25
    # 设置坐标轴标签和范围
    ax.set_xlabel('X', fontweight='bold', fontsize=14)
    ax.set_ylabel('Y', fontweight='bold', fontsize=14)
    ax.set_zlabel('Z', fontweight='bold', fontsize=14)
    ax.set_xlim([-20, 20])
    ax.set_ylim([-20, 20])
    ax.set_zlim([-10, 10])
    plt.title('Parametrinis blynas', fontweight='bold', fontsize=16)

    # 绘制三个正交坐标轴
    current_point = np.array([x[t-1], y[t-1], z[t-1]])
    ax.quiver(current_point[0], current_point[1], current_point[2], 
              emp * T[t, 0], emp * T[t, 1], emp * T[t, 2], color='r')  # 红色切向量T
    ax.quiver(current_point[0], current_point[1], current_point[2], 
              emp * B[t, 0], emp * B[t, 1], emp * B[t, 2], color='b')  # 蓝色副法向量B
    ax.quiver(current_point[0], current_point[1], current_point[2], 
              emp * N[t, 0], emp * N[t, 1], emp * N[t, 2], color='g')  # 绿色法向量N

    # 生成并绘制垂直于B的圆
    # 用T和N作为平面的正交基，将基础圆转换到当前平面
    rotation_matrix = np.array([T[t], N[t], B[t]]).T  # 正交旋转矩阵
    circle_3d = rotation_matrix @ circle_base
    # 平移到当前点
    circle_x = circle_3d[0] + current_point[0]
    circle_y = circle_3d[1] + current_point[1]
    circle_z = circle_3d[2] + current_point[2]
    # 绘制圆
    ax.plot(circle_x, circle_y, circle_z, color='cyan', linewidth=2)

    plt.pause(0.01)  # 控制动画播放速度
plt.show()

关键修改说明 #

• 提前生成圆的基础参数点circle_base，减少循环内重复计算，提升运行效率
• 利用正交标架T、N、B构建旋转矩阵，将平面圆转换到3D空间中垂直于B的平面
• 将生成的圆平移到当前轨迹点，确保圆随点移动且倾斜角度完全匹配B向量方向
• 增加导数归一化时的零向量判断，避免出现除以零的运行错误
• 调整循环起始值为1，防止t-1出现负索引越界问题

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