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

Question

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

Answer 1

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