To compute the derivatives, you need to use the same logic as in my previous answer:
```
import numpy as np
import matplotlib.pyplot as plt
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_values = np.linspace(0, 1, 50)
w_values = np.linspace(0, 1, 50)
u_grid, w_grid = np.meshgrid(u_values, w_values)
Q = np.zeros((u_grid.shape[0], u_grid.shape[1], 3))
Qu = np.zeros((u_grid.shape[0], u_grid.shape[1], 3))
Qw = np.zeros((u_grid.shape[0], u_grid.shape[1], 3))
for i in range(u_grid.shape[0]):
for j in range(u_grid.shape[1]):
u = u_grid[i, j]
w = w_grid[i, j]
U = np.array([u**3, u**2, u, 1])
W = np.array([[w**3], [w**2], [w], [1]])
U1 = np.array([3*u**2, 2*u, 1, 0])
W1 = np.array([[3*w**2], [2*w], [1], [0]])
Q[i, j, :] = np.array([U @ B_transformed[:, :, k] @ W for k in range(3)]).flatten()
Qu[i, j, :] = np.array([U1 @ B_transformed[:, :, k] @ W for k in range(3)]).flatten()
Qw[i, j, :] = np.array([U @ B_transformed[:, :, k] @ W1 for k in range(3)]).flatten()
print(" Q(0.5, 0.5):")
print(Q)
print("Qu(0.5, 0.5):")
print(Qu)
print("Qw(0.5, 0.5):")
print(Qw)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(Q[:, :, 0], Q[:, :, 1], Q[:, :, 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()
```
which returns:
```
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.]]]
Q(0.5, 0.5):
[[[-15. 0. 15. ]
[-14.3877551 0.29987505 15. ]
[-13.7755102 0.58725531 15. ]
...
[ 13.7755102 0.58725531 15. ]
[ 14.3877551 0.29987505 15. ]
[ 15. 0. 15. ]]
[[-15. 0.29987505 14.3877551 ]
[-14.3877551 0.58176509 14.3877551 ]
[-13.7755102 0.85190972 14.3877551 ]
...
[ 13.7755102 0.85190972 14.3877551 ]
[ 14.3877551 0.58176509 14.3877551 ]
[ 15. 0.29987505 14.3877551 ]]
[[-15. 0.58725531 13.7755102 ]
[-14.3877551 0.85190972 13.7755102 ]
[-13.7755102 1.10553686 13.7755102 ]
...
[ 13.7755102 1.10553686 13.7755102 ]
[ 14.3877551 0.85190972 13.7755102 ]
[ 15. 0.58725531 13.7755102 ]]
...
[[-15. 0.58725531 -13.7755102 ]
[-14.3877551 0.85190972 -13.7755102 ]
[-13.7755102 1.10553686 -13.7755102 ]
...
[ 13.7755102 1.10553686 -13.7755102 ]
[ 14.3877551 0.85190972 -13.7755102 ]
[ 15. 0.58725531 -13.7755102 ]]
[[-15. 0.29987505 -14.3877551 ]
[-14.3877551 0.58176509 -14.3877551 ]
[-13.7755102 0.85190972 -14.3877551 ]
...
[ 13.7755102 0.85190972 -14.3877551 ]
[ 14.3877551 0.58176509 -14.3877551 ]
[ 15. 0.29987505 -14.3877551 ]]
[[-15. 0. -15. ]
[-14.3877551 0.29987505 -15. ]
[-13.7755102 0.58725531 -15. ]
...
[ 13.7755102 0.58725531 -15. ]
[ 14.3877551 0.29987505 -15. ]
[ 15. 0. -15. ]]]
Qu(0.5, 0.5):
[[[ 30. 15. 0. ]
[ 30. 14.3877551 0. ]
[ 30. 13.7755102 0. ]
...
[ 30. -13.7755102 0. ]
[ 30. -14.3877551 0. ]
[ 30. -15. 0. ]]
[[ 30. 14.10037484 0. ]
[ 30. 13.52484934 0. ]
[ 30. 12.94932384 0. ]
...
[ 30. -12.94932384 0. ]
[ 30. -13.52484934 0. ]
[ 30. -14.10037484 0. ]]
[[ 30. 13.23823407 0. ]
[ 30. 12.69789798 0. ]
[ 30. 12.1575619 0. ]
...
[ 30. -12.1575619 0. ]
[ 30. -12.69789798 0. ]
[ 30. -13.23823407 0. ]]
...
[[ 30. 13.23823407 0. ]
[ 30. 12.69789798 0. ]
[ 30. 12.1575619 0. ]
...
[ 30. -12.1575619 0. ]
[ 30. -12.69789798 0. ]
[ 30. -13.23823407 0. ]]
[[ 30. 14.10037484 0. ]
[ 30. 13.52484934 0. ]
[ 30. 12.94932384 0. ]
...
[ 30. -12.94932384 0. ]
[ 30. -13.52484934 0. ]
[ 30. -14.10037484 0. ]]
[[ 30. 15. 0. ]
[ 30. 14.3877551 0. ]
[ 30. 13.7755102 0. ]
...
[ 30. -13.7755102 0. ]
[ 30. -14.3877551 0. ]
[ 30. -15. 0. ]]]
Qw(0.5, 0.5):
[[[ 0. 15. -30. ]
[ 0. 14.10037484 -30. ]
[ 0. 13.23823407 -30. ]
...
[ 0. 13.23823407 -30. ]
[ 0. 14.10037484 -30. ]
[ 0. 15. -30. ]]
[[ 0. 14.3877551 -30. ]
[ 0. 13.52484934 -30. ]
[ 0. 12.69789798 -30. ]
...
[ 0. 12.69789798 -30. ]
[ 0. 13.52484934 -30. ]
[ 0. 14.3877551 -30. ]]
[[ 0. 13.7755102 -30. ]
[ 0. 12.94932384 -30. ]
[ 0. 12.1575619 -30. ]
...
[ 0. 12.1575619 -30. ]
[ 0. 12.94932384 -30. ]
[ 0. 13.7755102 -30. ]]
...
[[ 0. -13.7755102 -30. ]
[ 0. -12.94932384 -30. ]
[ 0. -12.1575619 -30. ]
...
[ 0. -12.1575619 -30. ]
[ 0. -12.94932384 -30. ]
[ 0. -13.7755102 -30. ]]
[[ 0. -14.3877551 -30. ]
[ 0. -13.52484934 -30. ]
[ 0. -12.69789798 -30. ]
...
[ 0. -12.69789798 -30. ]
[ 0. -13.52484934 -30. ]
[ 0. -14.3877551 -30. ]]
[[ 0. -15. -30. ]
[ 0. -14.10037484 -30. ]
[ 0. -13.23823407 -30. ]
...
[ 0. -13.23823407 -30. ]
[ 0. -14.10037484 -30. ]
[ 0. -15. -30. ]]]
```
and the following plot
[![enter image description here][1]][1]
[1]: https://i.sstatic.net/TSqA8RJj.png
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
[![enter image description here][1]][1]
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
[![enter image description here][2]][2]
[1]: https://i.sstatic.net/V0zMfjZt.png
[2]: https://i.sstatic.net/E4dConQZ.png
