Question: 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

Answer: 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 your particular problem, you can do this: 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