Bezier曲面矩阵构造问题：矩阵乘法代码错误排查

问题背景 #

需要计算Bezier曲面的矩阵变换公式 N @ B @ N^T，其中控制点矩阵B的正确定义为：

import numpy as np

# 书籍给出的控制点矩阵B
B = np.array([
    [[-15, 0, 15], [-15, 5, 5], [-15, 5, -5], [-15, 0, -15]], 
    [[-5, 5, 15], [-5, 5, 5], [-5, 5, -5], [-5, 5, -15]], 
    [[5, 5, 15], [5, 5, 5], [5, 5, -5], [5, 5, -15]], 
    [[15, 0, 15], [15, 5, 5], [15, 5, -5], [15, 0, -15]]
])

变换矩阵N为：

N = np.array([
    [-1, 3, -3, 1],
    [3, -6, 3, 0],
    [-3, 3, 0, 0],
    [1, 0, 0, 0]
])
Nt = N.T  # 正确的转置，无需手动重复定义

自行编写的代码逻辑错误，得到的结果与书籍给出的正确结果不符，同时后续计算曲面点Q(0.5, 0.5)也无法得到预期值[0, 4.6875, 0]。

错误代码分析 #

原代码中的三重循环完全误解了矩阵乘法的核心逻辑：

B_transformed = np.zeros_like(B)
for i in range(B.shape[0]):
    for j in range(B.shape[1]): 
        for k in range(3): 
            B_transformed[i, j, k] = B[i, j, k] * N[j, k] * Nt[j, k]

• 矩阵乘法的核心是行与列对应元素相乘后求和，而非单个元素的直接相乘
• 控制点矩阵B是4×4×3的张量，每个维度（x/y/z坐标）需要单独作为4×4矩阵执行N @ B_channel @ Nt的运算

正确实现代码 #

步骤1：计算N @ B @ N^T #

# 初始化结果张量
NBNt = np.zeros_like(B)

# 对每个坐标通道（x/y/z）分别执行矩阵乘法
for channel in range(3):
    # 提取当前通道的4×4矩阵
    B_channel = B[:, :, channel]
    # 执行矩阵乘法：N × B_channel × N^T
    NBNt[:, :, channel] = N @ B_channel @ Nt

# 输出结果，与书籍给出的正确值一致
print(NBNt)

运行后得到的结果：

[[[0 0 0]
  [0 0 0]
  [0 0 0]
  [0 0 0]]

 [[0 0 0]
  [0 -45 0]
  [0 45 0]
  [0 -15 0]]

 [[0 0 0]
  [0 45 0]
  [0 -45 0]
  [30 15 0]]

 [[0 0 0]
  [0 -15 0]
  [0 15 -30]
  [-15 0 15]]]

步骤2：计算曲面点Q(0.5, 0.5) #

根据公式U @ NBNt @ U^T（其中U = [0.125, 0.25, 0.5, 1.0]），代码实现如下：

U = np.array([0.125, 0.25, 0.5, 1.0])
# 计算U × NBNt × U^T，得到3维坐标点
Q = U @ NBNt @ U.reshape(-1, 1)
# 转换为一维数组输出
print(Q.flatten())

运行结果为：[0. 4.6875 0. ]，与预期一致。

关键逻辑总结 #

• 张量形式的控制点矩阵需要按坐标通道拆分，每个通道单独执行矩阵乘法
• 矩阵乘法必须遵循行乘列求和的规则，numpy的@运算符会自动处理该逻辑
• 转置矩阵直接使用.T属性获取，避免手动输入导致的错误

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