Matplotlib 3D多边形旋转后发生扭曲的问题求助

我现在在做一个3D物体旋转的功能：通过指定物体的一个面，让这个面的法向量对齐到目标向量。如果目标向量是平行于x、y或者z轴的，旋转出来的效果完全正常，但只要换成任意非轴向的向量，整个物体就会变得扭曲变形。

比如我代码里第104行的那个目标向量，运行后就能看到扭曲的问题。我猜测是不是需要额外加个偏移或者修正，但完全摸不着头绪。有没有大佬能帮我找找问题出在哪呀？

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import numpy as np

# Create a new figure with 3D projection
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')


def rotate_shape_to_normal(vertices, face_indices, target_normal):
    """
    Rotate vertices so that the face specified by face_indices becomes normal to target_normal.

    :param vertices: List or numpy array of 3D points defining the shape
    :param face_indices: Indices of the vertices that form the face to align
    :param target_normal: The vector to which the face should be normal
    :return: Rotated vertices
    """
    vertices = np.array(vertices)  # Convert to numpy array if it isn't already
    face_vertices = vertices[face_indices]

    # Calculate the normal vector of the face
    v1 = face_vertices[1] - face_vertices[0]
    v2 = face_vertices[2] - face_vertices[0]
    face_normal = np.cross(v1, v2)
    print(face_normal)
    ax.quiver(*face_vertices[0], *face_normal, color='r', length=2)
    ax.quiver(*face_vertices[0], *target_normal, color='b', length=2)
    face_normal = face_normal / np.linalg.norm(face_normal)  # Normalize

    # Normalize target_normal for consistency
    target_normal = target_normal / np.linalg.norm(target_normal)

    # Compute the rotation axis (perpendicular to both vectors)
    rotation_axis = np.cross(face_normal, target_normal)
    rotation_axis = rotation_axis / np.linalg.norm(rotation_axis)# if np.linalg.norm(rotation_axis) > 0 else np.array([1, 0, 0])  # Default to x-axis if parallel
    print("rotation_axis", rotation_axis)
    # Compute the angle between the vectors
    cos_theta = np.dot(face_normal, target_normal)
    print("cos_theta", np.degrees(cos_theta))
    theta = np.arccos(np.clip(cos_theta, -1.0, 1.0))  # Clip to avoid floating-point issues
    print("theta", np.degrees(theta))

    # Check if vectors are already aligned or opposite
    if np.isclose(theta, 0) or np.isclose(theta, np.pi):
        return vertices  # No rotation needed

    # Rodrigues' rotation formula
    K = np.array([
        [0, -rotation_axis[2], rotation_axis[1]],
        [rotation_axis[2], 0, -rotation_axis[0]],
        [-rotation_axis[1], rotation_axis[0], 0]
    ])
    rotation_matrix = np.sin(theta) * K + (1 - np.cos(theta)) * np.outer(rotation_axis, rotation_axis)


    # Apply the rotation to all vertices
    return np.dot(vertices, rotation_matrix.T)


# Define the vertices for a simple satellite model:
vertices = np.array([
    # Main body (cube)
    (0, 0, 0),  # 0 - Base, front-left
    (1, 0, 0),  # 1 - Base, front-right
    (1, 1, 0),  # 2 - Base, back-right
    (0, 1, 0),  # 3 - Base, back-left

    (0, 0, 1),  # 4 - Top, front-left
    (1, 0, 1),  # 5 - Top, front-right
    (1, 1, 1),  # 6 - Top, back-right
    (0, 1, 1),  # 7 - Top, back-left

    # Solar panels
    # Left panel
    (0.5, .25, -1.5), # 8
    (0.5, .75, -1.5), # 9
    (0.5, .75, 0), # 10
    (0.5, .25, 0), # 11

    # Right panel
    (0.5, .25, 2.5),  # 12
    (0.5, .75, 2.5),  # 13
    (0.5, .75, 1),  # 14
    (0.5, .25, 1)   # 15

])



# Translate spacecraft to origin
vertices = [(x - .5, y - .5, z - .5) for x, y, z in vertices]

# Original points
ax.scatter(vertices[5][0], vertices[5][1],vertices[5][2], color='green', s=25)
ax.scatter(vertices[1][0], vertices[1][1],vertices[1][2], color='green', s=25)
ax.scatter(vertices[4][0], vertices[4][1],vertices[4][2], color='green', s=25)


# Rotate spacecraft to point at specified vector
face_to_align = [5, 1, 4] # Indexes of points that define the face I want to be normal to the specified vector (-Y)
#vector = np.array([0, 0, 1])  # This rotation doesn't skew anything
vector = np.array([0.56167836, 0.76075023, 0.32523301]) # This rotation skews everything

vertices = rotate_shape_to_normal(vertices, face_to_align, vector)

# Rotated points
ax.scatter(vertices[5][0], vertices[5][1],vertices[5][2], color='blue', s=25)
ax.scatter(vertices[1][0], vertices[1][1],vertices[1][2], color='blue', s=25)
ax.scatter(vertices[4][0], vertices[4][1],vertices[4][2], color='blue', s=25)


# Faces for the cube
cube_faces = [
    [vertices[0], vertices[1], vertices[5], vertices[4]],  # (-Y)
    [vertices[1], vertices[2], vertices[6], vertices[5]],  # (+X)
    [vertices[2], vertices[3], vertices[7], vertices[6]],  # (+Y)
    [vertices[3], vertices[0], vertices[4], vertices[7]],  # (-X)
    [vertices[0], vertices[1], vertices[2], vertices[3]],  # (-Z)
    [vertices[4], vertices[5], vertices[6], vertices[7]],  # (+Z)
]

# Solar panel faces - note these are just rectangles
left_panel = [vertices[8], vertices[9], vertices[10], vertices[11]]
right_panel = [vertices[12], vertices[13], vertices[14], vertices[15]]

# Combine all faces
faces = cube_faces + [left_panel, right_panel]


# Create Poly3DCollection
poly3d = Poly3DCollection(faces, alpha=0.7)
poly3d.set_edgecolor('k')

# Set face color for different parts
poly3d.set_facecolor([[1, 0, 0], [.7, .7, .7], [.7, .7, .7], [.7, .7, .7], [.7, .7, .7], [.7, .7, .7], [0, 0.5, 1], [0, 0.5, 1]])
ax.add_collection3d(poly3d)

# Set the aspect ratio to ensure it looks like a cube
ax.set_box_aspect((1, 1, 1))

# Remove axes for a cleaner look
#ax.set_axis_off()

# Set limits to see everything
ax.set_xlim(-1, 2)
ax.set_ylim(-1, 2)
ax.set_zlim(-1, 2)

# Origin
ax.scatter(0, 0, 0, color='red', s=25)


ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')

# Adjust view for better visualization
ax.view_init(elev=20., azim=-45)

plt.show()

备注：内容来源于stack exchange，提问作者dreed75