I am trying to rewrite slow method("bezier_surf_eval") for bezier surface with quick one(in matrix form "bezier_surface_M").
Using this formula:
Were [N] and [M] an bezier basis matrix (4x2 degree, because of 5x3 control points), [U], [W] is surface uv params, [B] - control points
What i do wrong?
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm
from matplotlib.figure import Figure
from mpl_toolkits.mplot3d.axes3d import Axes3D
from scipy.special import comb
def bern(i, d, t):
return comb(d, i) * (t ** (d - i)) * (1 - t) ** i
def bezier_matrix(d):
return np.array([[(-1) ** (i - j) * comb(j, i) * comb(d, j) for i in range(d 1)] for j in range(d 1)], int)
BM = [bezier_matrix(i) for i in range(16)]
poles5x3 = np.array([
[[1.8, -0.3, 0.], [1.8, 0.13, 0.1], [1.8, 0.5, 0.]],
[[2., -0.3, 0.06], [2.1, 0.1, 0.1], [2.1, 0.5, 0.1]],
[[2.3, -0.3, 0.1], [2.3, 0.13, 0.2], [2.3, 0.5, 0.1]],
[[2.4, -0.3, 0.1], [2.5, 0.1, 0.15], [2.5, 0.5, 0.1]],
[[2.6, -0.3, 0.], [2.6, 0.1, 0.1], [2.5, 0.5, 0.]]])
def bezier_surf_eval(poles, u, v, deg_u, deg_v):
count_u, count_v = deg_u 1, deg_v 1
co = [0, 0, 0]
for i in range(count_u):
for j in range(count_v):
bu = bern(i, deg_u, u)
bv = bern(j, deg_v, v)
pole_co = poles[i][j]
m = bu * bv
co[0] = pole_co[0] * m
co[1] = pole_co[1] * m
co[2] = pole_co[2] * m
return co
def bezier_surface_slow(poles, resol=(16, 16)):
params_u = np.linspace(0, 1, resol[0])
params_v = np.linspace(0, 1, resol[1])
count_u, count_v, _ = poles.shape
du, dv = count_u - 1, count_v - 1
cps = poles.tolist()
coords = np.empty((*resol, 3), float)
for vi, v in enumerate(params_v):
for ui, u in enumerate(params_u):
coords[ui, vi] = bezier_surf_eval(cps, u, v, deg_u=du, deg_v=dv)
return np.array(coords, np.float32)
def bezier_surface_M(cps: np.ndarray, resol=(16, 16)):
u = np.linspace(0, 1, resol[0])
v = np.linspace(0, 1, resol[1])
count_u, count_v, _ = cps.shape
deg_u, deg_v = count_u - 1, count_v - 1
u_vec = np.array([u ** i for i in range(count_u)])
v_vec = np.array([v ** i for i in range(count_v)])
BM_u, BM_v = BM[deg_u], BM[deg_v]
return u_vec.T.dot(BM_u).dot(cps).dot(BM_v.T).dot(v_vec)
fig: Figure = plt.figure(figsize=(7, 7))
ax: Axes3D = fig.add_subplot(111, projection='3d')
# ax.scatter(px, py, pz)
# --------------
resol = 16, 16
co = bezier_surface_slow(poles5x3, resol=resol)
x, y, z = co.T
ax.plot_surface(x, y, z, cmap=cm.gray, linewidth=1, antialiased=False)
# --------------
resol = 16, 16
co = bezier_surface_M(poles5x3, resol=resol)
x, y, z = co.T
ax.plot_surface(x, y, z, cmap=cm.gray, linewidth=1, antialiased=False)
plt.show()
CodePudding user response:
Looks like i solved it. Problem was in bad order for matrices (in method "bezier_surface_M") for dot product. Multiplication with per-axis matrices with correct orders do the job.
Replaced:
u_vec.T.dot(BM_u).dot(cps).dot(BM_v.T).dot(v_vec)
With:
cps_x = cps[:, :, 0]
cps_y = cps[:, :, 1]
cps_z = cps[:, :, 2]
m1 = u_vec.T.dot(BM_u)
m2 = BM_v.T.dot(v_vec)
x = m1.dot(cps_x).dot(m2)
y = m1.dot(cps_y).dot(m2)
z = m1.dot(cps_z).dot(m2)
Maybe it can be optimized with np.tensorproduct or something like that, but i`m not as good with numpy and matrices math.
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm
from matplotlib.figure import Figure
from mpl_toolkits.mplot3d.axes3d import Axes3D
from scipy.special import comb
def bezier_matrix(d):
return np.array([[(-1) ** (i - j) * comb(j, i) * comb(d, j) for i in range(d 1)] for j in range(d 1)], int)
BM = [bezier_matrix(i) for i in range(16)]
poles5x3 = np.array([
[[1.8, -0.3, 0.], [1.8, 0.13, 0.1], [1.8, 0.5, 0.]],
[[2., -0.3, 0.06], [2.1, 0.1, 0.1], [2.1, 0.5, 0.1]],
[[2.3, -0.3, 0.1], [2.3, 0.13, 0.2], [2.3, 0.5, 0.1]],
[[2.4, -0.3, 0.1], [2.5, 0.1, 0.15], [2.5, 0.5, 0.1]],
[[2.6, -0.3, 0.], [2.6, 0.1, 0.1], [2.5, 0.5, 0.]]])
def bezier_surface_M(cps: np.ndarray, resol=(16, 16)):
u, v = np.linspace(0, 1, resol[0]), np.linspace(0, 1, resol[1])
count_u, count_v, _ = cps.shape
deg_u, deg_v = count_u - 1, count_v - 1
u_vec = np.array([u ** i for i in range(count_u)])
v_vec = np.array([v ** i for i in range(count_v)])
BM_u, BM_v = BM[deg_u], BM[deg_v]
cps_x = cps[:, :, 0]
cps_y = cps[:, :, 1]
cps_z = cps[:, :, 2]
m1 = u_vec.T.dot(BM_u)
m2 = BM_v.T.dot(v_vec)
x = m1.dot(cps_x).dot(m2)
y = m1.dot(cps_y).dot(m2)
z = m1.dot(cps_z).dot(m2)
return x, y, z
fig: Figure = plt.figure(figsize=(7, 7))
ax: Axes3D = fig.add_subplot(111, projection='3d')
resol = 32, 32
# --------------
x, y, z = bezier_surface_M(poles5x3, resol=resol)
ax.plot_surface(x, y, z, cmap=cm.gray, linewidth=1, antialiased=False)
plt.show()