矩阵形式的 Bezier MxN 曲面（使用 Python 和 numpy）

Question

我正在尝试用快速方法（矩阵形式“bezier_surface_M”）为贝塞尔曲面重写慢速方法（“bezier_surf_eval”）。

使用这个公式：

[N] 和 [M] 是贝塞尔曲线基础矩阵（4x2 度，因为有 5x3 个控制点）， [U]、[W] 是表面 uv 参数，[B] - 控制点

我做错了什么？

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()

Answer 1

看来我解决了。问题是点积矩阵的顺序不正确（在方法“bezier_surface_M”中）。与具有正确顺序的 per-axis 矩阵相乘即可。

已替换：

u_vec.T.dot(BM_u).dot(cps).dot(BM_v.T).dot(v_vec)

与：

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)

也许它可以用 np.tensorproduct 或类似的东西优化，但我对 numpy 和矩阵数学不太好。

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()

矩阵形式的 Bezier MxN 曲面（使用 Python 和 numpy）

Bezier MxN Surface in matrix form( with Python and numpy)

python

numpy

bezier

surface