Матричная форма поверхности Безье

У меня возникла проблема с построением поверхности Безье по примеру из книги с использованием математических формул в матричной форме. Особенно при умножении матриц.

Я пытаюсь использовать эту формулу У меня есть матрица контрольных точек

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][B][N]^t

И я попробовал умножить матрицу на эти два, но у итоговой матрицы получаю совсем другие значения, я понимаю, что скорее всего проблема в коде

"

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]]]

Правильный ответ из книги:

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

Далее также будет производиться умножение матриц, поэтому мне важно понять, что я делаю не так.

Q(0,5, 0,5) =

[0.125 0.25  0.5   1.   ] * [N][B][N]^t * [[0.125]
                                           [0.25 ]
                                           [0.5  ]
                                           [1.   ]]

Это расчет точки на поверхности при w = 0,5 и u = 0,5.

И ответ должен быть
[0, 4,6875, 0]

Я использую блокнот Jupyter

Добро пожаловать в ТАК! Действительно хороший первый вопрос! Приятного кодирования!

Serge de Gosson de Varennes 06.06.2024 17:11
Почему в Python есть оператор "pass"?
Почему в Python есть оператор "pass"?
Оператор pass в Python - это простая концепция, которую могут быстро освоить даже новички без опыта программирования.
Некоторые методы, о которых вы не знали, что они существуют в Python
Некоторые методы, о которых вы не знали, что они существуют в Python
Python - самый известный и самый простой в изучении язык в наши дни. Имея широкий спектр применения в области машинного обучения, Data Science,...
Основы Python Часть I
Основы Python Часть I
Вы когда-нибудь задумывались, почему в программах на Python вы видите приведенный ниже код?
LeetCode - 1579. Удаление максимального числа ребер для сохранения полной проходимости графа
LeetCode - 1579. Удаление максимального числа ребер для сохранения полной проходимости графа
Алиса и Боб имеют неориентированный граф из n узлов и трех типов ребер:
Оптимизация кода с помощью тернарного оператора Python
Оптимизация кода с помощью тернарного оператора Python
И последнее, что мы хотели бы показать вам, прежде чем двигаться дальше, это
Советы по эффективной веб-разработке с помощью Python
Советы по эффективной веб-разработке с помощью Python
Как веб-разработчик, Python может стать мощным инструментом для создания эффективных и масштабируемых веб-приложений.
2
1
55
1
Перейти к ответу Данный вопрос помечен как решенный

Ответы 1

Ответ принят как подходящий

Обычно поверхность Безье строится таким образом (как вопрос опубликован в 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()

которые возвращаются

Теперь для вашей конкретной проблемы вы можете сделать это:

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)

что дает вам

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.    ]

и если вы также хотите построить это, вы можете адаптировать мой верхний код к этому:

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

дарю тебе снова

Другие вопросы по теме