Bezier 曲面矩阵形式

Question

我在按照书中的示例使用矩阵形式的数学公式构造贝塞尔曲面时遇到问题。 尤其是矩阵相乘时。

我正在尝试使用这个公式 我有一个控制点矩阵

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][B][N]^t

我尝试将矩阵乘以这两个，但我得到的最终矩阵完全不同的值，我知道问题很可能出在代码中

”

B = np.array([
    [[-15, 0, 15], [-5, 5, 15], [5, 5, 15], [15, 0, 15]], 
    [[-15, 5, 5], [-5, 5, 5], [5, 5, 5], [15, 5, 5]], 
    [[-15, 5, -5], [-5, 5, -5], [5, 5, -5], [15, 5, -5]], 
    [[-15, 0, -15], [-5, 5, -15], [5, 5, -15], [15, 0, -15]]
])

N = np.array([[-1, 3, -3, 1],
               [3, -6, 3, 0],
               [-3, 3, 0, 0],
               [1, 0, 0, 0]
              ])

Nt = np.array([[-1, 3, -3, 1],
               [3, -6, 3, 0],
               [-3, 3, 0, 0],
               [1, 0, 0, 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]

”

 [[[ -15    0  135]
  [ -45  180  135]
  [  45   45    0]
  [  15    0    0]]

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

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

 [[ -15    0 -135]
  [ -45  180 -135]
  [  45   45    0]
  [  15    0    0]]]

书中正确答案是

NBNt = np.array([
    [[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]]
])

接下来，还将执行矩阵乘法，所以了解我做错了什么对我来说很重要

Q(0.5, 0.5) =

[0.125 0.25  0.5   1.   ] * [N][B][N]^t * [[0.125]
                                           [0.25 ]
                                           [0.5  ]
                                           [1.   ]]

这是表面上 w = 0.5 且 u = 0.5 处的点的计算

答案应该是
[0, 4.6875, 0]

我使用 Jupyter Notebook

Answer 1

通常，贝塞尔曲面是这样绘制的（正如问题发布在

matplotlib

中）。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.special import comb

def bernstein_poly(i, n, t):
    return comb(n, i) * (t**i) * ((1 - t)**(n - i))

def bernstein_matrix(n, t):
    return np.array([bernstein_poly(i, n, t) for i in range(n + 1)])

P = 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, m = P.shape[0] - 1, P.shape[1] - 1

u = np.linspace(0, 1, 50)
v = np.linspace(0, 1, 50)
U, V = np.meshgrid(u, v)

surface_points = np.zeros((U.shape[0], U.shape[1], 3))
for i in range(U.shape[0]):
    for j in range(U.shape[1]):
        Bu = bernstein_matrix(n, U[i, j])
        Bv = bernstein_matrix(m, V[i, j])
        surface_points[i, j] = np.tensordot(np.tensordot(Bu, P, axes=(0, 0)), Bv, axes=(0, 0))

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(surface_points[:,:,0], surface_points[:,:,1], surface_points[:,:,2], rstride=1, cstride=1, color='b', alpha=0.6, edgecolor='w')
ax.scatter(P[:,:,0], P[:,:,1], P[:,:,2], color='r', s=50)

plt.show()

返回

现在，对于您的特定问题，您可以这样做：

import numpy as np

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 = np.array([[-1, 3, -3, 1],
              [3, -6, 3, 0],
              [-3, 3, 0, 0],
              [1, 0, 0, 0]])

Nt = N.T

B_transformed = np.zeros((4, 4, 3))

for i in range(3):  
    B_transformed[:, :, i] = N @ B[:, :, i] @ Nt

print("Transformed control points matrix B_transformed:")
print(B_transformed)

u = 0.5
w = 0.5

U = np.array([u**3, u**2, u, 1])
W = np.array([w**3, w**2, w, 1])

Q = np.array([U @ B_transformed[:, :, i] @ W for i in range(3)])

print("Point on the Bézier surface Q(0.5, 0.5):")
print(Q)

这给了你

Transformed control points matrix B_transformed:
[[[  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.]]]
Point on the Bézier surface Q(0.5, 0.5):
[0.     4.6875 0.    ]

如果你也想绘制它，你可以调整我的顶部代码：

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.special import comb

def bernstein_poly(i, n, t):
    return comb(n, i) * (t**i) * ((1 - t)**(n - i))

def bernstein_matrix(n, t):
    return np.array([bernstein_poly(i, n, t) for i in range(n + 1)])

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 = np.array([[-1, 3, -3, 1],
              [3, -6, 3, 0],
              [-3, 3, 0, 0],
              [1, 0, 0, 0]])

Nt = N.T

B_transformed = np.zeros((4, 4, 3))

for i in range(3): 
    B_transformed[:, :, i] = N @ B[:, :, i] @ Nt

print("Transformed control points matrix B_transformed:")
print(B_transformed)

u = np.linspace(0, 1, 50)
w = np.linspace(0, 1, 50)
U, W = np.meshgrid(u, w)

surface_points = np.zeros((U.shape[0], U.shape[1], 3))
for i in range(U.shape[0]):
    for j in range(U.shape[1]):
        U_vec = np.array([U[i, j]**3, U[i, j]**2, U[i, j], 1])
        W_vec = np.array([W[i, j]**3, W[i, j]**2, W[i, j], 1])
        surface_points[i, j] = np.array([U_vec @ B_transformed[:, :, k] @ W_vec for k in range(3)])

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(surface_points[:,:,0], surface_points[:,:,1], surface_points[:,:,2], rstride=1, cstride=1, color='b', alpha=0.6, edgecolor='w')
ax.scatter(B[:,:,0], B[:,:,1], B[:,:,2], color='r', s=50)

plt.show()

再次给你