I have a problem with constructing a Bezier surface following an example from a book, using mathematical formulas in matrix form.
Especially when multiplying matrices.
I’m trying to use this formula
I have a matrix of control points
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]]
])
And we have to multiply it by matrices
and get [N][B][N]^t
And I tried to multiply the matrix by these two, but I get completely different values ​​for the final matrix, I understand that most likely the problem is in the code
"
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]]]
Correct answer from book is
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]]
])
Next, matrix multiplication will also be performed, so it’s important for me to understand what I’m doing wrong
Q(0.5, 0.5) =
[0.125 0.25 0.5 1. ] * [N][B][N]^t * [[0.125]
[0.25 ]
[0.5 ]
[1. ]]
This is the calculation of a point on a surface at w = 0.5 and u = 0.5
And the answer should be
[0, 4.6875, 0]
I use Jupyter Notebook
>Solution :
Generally, Bezier surface are plotted this way (as the question is posted in 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()
which return
Now, for you particular problem, you can dothis:
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)
which gives you
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. ]
and if you also want to plot it, you can adapt my top code to this:
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()
giving you again
