De certaines circonstances, j'étudiais la division de la courbe de Bézier. Je n'ai pas trouvé beaucoup de sites expliquant comment diviser la ** courbe de Bézier d'ordre $ n $ **.
Eh bien, puisque la courbe de Bézier n'est utilisée que jusqu'au 4ème ordre, peu importe si vous ne connaissez pas l'ordre $ n $. Je voulais vraiment le mettre en œuvre proprement, alors je l'ai étudié, donc je vais le laisser comme connaissance. J'ai posté un exemple d'implémentation dans l'article, mais je l'ai également posté sur GitHub (voir line_module.py).
Il s'agit d'une méthode appelée algorithme de De Casteljau. Cette méthode utilise la récurrence pour trouver les points de contrôle de la nouvelle courbe de Bézier après la division. Cette méthode est ** de conception simple et facile à mettre en œuvre **, mais présente l'inconvénient d'être difficile à utiliser dans les systèmes embarqués avec une petite taille de pile **. Eh bien, la courbe de Bézier n'est que jusqu'au 4ème ordre (omis)
À titre d'exemple, divisons la courbe de Bézier ci-dessous. Soit $ t $ la variable médiatrice de la courbe de Bézier au point de division, et soit $ P_i $ le point de contrôle initial de la courbe de Bézier.
** 1. Trouvez le point qui devient $ t: 1-t $ sur la ligne reliant les points de contrôle de la courbe de Bézier (point $ Q $). ** **
** 2. Pour la ligne reliant les points obtenus, trouvez le point qui est $ t: 1-t $ (point $ R $). ** **
** 3. Répétez jusqu'à ce qu'il y ait un point (point $ S $). ** **
** 4. À tous les points, y compris les points de contrôle initiaux, l'ensemble des points avec le plus petit indice et l'ensemble des points avec le plus grand indice sont les points de contrôle après la division. ** **
Pour la dernière partie (indice minimum), il sera plus facile à comprendre si vous écrivez l'arborescence suivante. L'arborescence ci-dessous montre à partir de quel segment de ligne (de quels deux points) un nouveau point a été généré.
Les points les plus extérieurs de cette arborescence, à savoir les points $ \ {P_0, Q_0, R_0, S_0 \} $ et les points $ \ {S_0, R_1, Q_2, P_3 \} $, deviennent les nouvelles courbes de Bézier. Je vais. Vous pouvez voir que les points de début et de fin de la courbe de Bézier d'origine sont communs et que le point (point de division $ S_0 $) dans la variable médiatrice $ t $ est le point de départ (point final) de la nouvelle courbe de Bézier.
Un exemple d'implémentation d'un programme fractionné utilisant la récurrence est présenté ci-dessous.
import numpy as np
from scipy.special import comb
import matplotlib.pyplot as plt
class Bezier:
def __init__(self, points):
points = [np.array(e).flatten() for e in points]
self._n = len(points)
self._points = points
@property
def dims(self):
return self._n - 1
@property
def points(self):
return self._points
@property
def points4matplot(self):
xs, ys = zip(*self.points)
return xs, ys
def _bernstein(self, n, i, t):
return comb(n, i) * t**i * (1 - t)**(n-i)
def bezier_point(self, t):
p = np.zeros(2)
for i, f in enumerate(self.points):
p += self._bernstein(self.dims, i, t) * f
return np.array(p).flatten()
def _de_casteljau_algorithm(self, points, t):
prev_p = None
q = []
for p in points:
if not prev_p is None:
## Split line into t:(1-t)
q.append(np.array((1-t) * prev_p + t * p).flatten())
prev_p = p
if len(q) == 1:
return [q]
return [q] + self._de_casteljau_algorithm(q, t)
def split_with_recursion(self, t):
ret = [self.points] + self._de_casteljau_algorithm(self.points, t)
lp = []
rp = []
for lst in ret:
## Fetch min and max point
lp.append(lst[0])
rp.append(lst[-1])
return Bezier(lp), Bezier(rp)
def plot(self, ax=None, with_ctrl_pt=True, bcolor="black", ccolor="gray", resolution=100):
if ax is None:
_, ax = plt.subplots()
prev_point = None
for t in np.linspace(0, 1, resolution):
bp = self.bezier_point(t)
if prev_point is None:
prev_point = (bp[0], bp[1])
xs, ys = zip(*(prev_point, (bp[0], bp[1])))
ax.plot(xs, ys, '-', color=bcolor)
prev_point = (bp[0], bp[1])
if with_ctrl_pt:
xs, ys = self.points4matplot
ax.plot(xs, ys, 'x--', lw=2, color=ccolor, ms=10)
return ax
def main():
bezier = Bezier([(0.3, 1), (0.2, 3), (0.4, 4), (0.5, 0)])
ax = bezier.plot()
left_bezier, right_bezier = bezier.split_with_recursion(0.3)
left_bezier.plot(ax, bcolor = "red")
right_bezier.plot(ax, bcolor = "blue")
plt.grid()
plt.show()
if __name__ == "__main__":
main()
Il s'agit d'une méthode pour trouver un nouveau point de contrôle à l'aide d'une matrice. Cette méthode ** n'a pas d'importance si la taille de la pile est petite **, mais ** le calcul est un peu compliqué **. Eh bien, la courbe de Bézier est au plus du 4e ordre (omis)
Comme étape préliminaire pour expliquer l'algorithme, j'expliquerai la représentation matricielle de la courbe de Bézier. La courbe de Bézier peut être exprimée par l'équation suivante.
\begin{align}
F(t) &= \sum_{i=0}^{n} P_i~B_{i}^{n}(t)
\end{align}
$ P_i $ représente le point de contrôle et $ B_ {i} ^ {n} (t) $ représente la fonction de base de Bernstein. $ n $ représente l'ordre de la courbe de Bézier plus $ 1 $ (le nombre de $ = $ points de contrôle). Cela peut être exprimé sous la forme d'un produit de matrices comme suit.
\begin{align}
F(t) &= \sum_{i=0}^{n} P_i~B_{i}^{n}(t) \\
&= \left(
\begin{array}{ccc}
B_{0}^{n} & B_{1}^{n} & \dots & B_{n}^{n}
\end{array}
\right)
\left(
\begin{array}{ccc}
P_0 \\
P_1 \\
\vdots \\
P_n
\end{array}
\right)
\end{align}
De plus, la formule de la courbe de Bézier peut être exprimée comme suit en l'organisant pour $ t $. Ici, $ a_n $ représente un coefficient approprié.
\begin{align}
F(t) &= \sum_{i=0}^{n} P_i~B_{i}^{n}(t) \\
&= \left(
\begin{array}{ccc}
1 & t & t^2 & \dots & t^n
\end{array}
\right)
\left(
\begin{array}{ccc}
a_0 \\
a_1 \\
\vdots \\
a_n
\end{array}
\right)
\end{align}
De ce qui précède, l'équation suivante peut être dérivée.
\begin{align}
F(t) &=
\left(
\begin{array}{ccc}
B_{0}^{n} & B_{1}^{n} & \dots & B_{n}^{n}
\end{array}
\right)
\left(
\begin{array}{ccc}
P_0 \\
P_1 \\
\vdots \\
P_n
\end{array}
\right)
=
\left(
\begin{array}{ccc}
1 & t & t^2 & \dots & t^n
\end{array}
\right)
\left(
\begin{array}{ccc}
a_0 \\
a_1 \\
\vdots \\
a_n
\end{array}
\right) \\ \\
F(t) &= BP = XA
\end{align}
C'est assez propre. Ici, nous savons que $ A $ peut être exprimé comme suit par une matrice triangulaire inférieure $ U_t $ et un point de contrôle $ P $.
\begin{align}
A &= U_tP \\ \\
&= \left(
\begin{array}{ccc}
{n \choose 0} {n \choose 0} (-1)^0 & 0 & \cdots & 0 \\
{n \choose 0} {n \choose 1} (-1)^1 & {n \choose 1} {n-1 \choose 0} (-1)^0 & \cdots & 0 \\
\vdots & \vdots & \cdots & \vdots \\
{n \choose 0} {n \choose n} (-1)^n & {n \choose 1} {n-1 \choose n-1} (-1)^{n-1} & \cdots & {n \choose n} {n-n \choose n-n} (-1)^{n-n} \\
\end{array}
\right)
\left(
\begin{array}{ccc}
P_0 \\
P_1 \\
\vdots \\
P_n
\end{array}
\right)
\end{align}
Par conséquent, finalement, l'équation de la courbe de Bézier peut être exprimée en utilisant une matrice comme suit.
F(t) = BP = XU_tP
L'algorithme de division de la courbe de Bézier à l'aide d'une matrice est le suivant.
** 1. Organisez la courbe de Bézier pour $ t $ et exprimez-la sous la forme d'une matrice. ** **
\begin{align}
F(t) &= XU_tP \\
&= \left(
\begin{array}{ccc}
1 & t & t^2 & \dots & t^n
\end{array}
\right)U_tP \\
\end{align}
** 2. Séparez le point de division $ z ~ (z \ in t) $. ** **
\begin{align}
F(t) &= \left(
\begin{array}{ccc}
1 & t & t^2 & \dots & t^n
\end{array}
\right)U_tP \\ \\
&= \left(
\begin{array}{ccc}
1 & (z\cdot t) & (z\cdot t)^2 & \dots & (z\cdot t)^n
\end{array}
\right)U_tP \\ \\
&= \left(
\begin{array}{ccc}
1 & t & t^2 & \dots & t^n
\end{array}
\right)
\left(
\begin{array}{ccc}
1 & & & & 0 \\
& z & & & \\
& & z^2 & & \\
& & & \ddots & \\
0 & & & & z^n
\end{array}
\right)U_tP
\end{align}
** 3. Transformez la formule comme suit pour trouver $ Q $. ** **
\begin{align}
F(t) &= \left(
\begin{array}{ccc}
1 & t & t^2 & \dots & t^n
\end{array}
\right)
\left(
\begin{array}{ccc}
1 & & & & 0 \\
& z & & & \\
& & z^2 & & \\
& & & \ddots & \\
0 & & & & z^n
\end{array}
\right)U_tP \\
&= X ~ Z U_t ~ P \\ \\
&= X ~ U_t U_t^{-1} Z U_t ~ P \\ \\
&= X ~ U_t Q ~P
\end{align}
Dans le calcul ci-dessus, la position de $ U_t $ est déplacée. Puisque $ U_t U_t ^ {-1} $ est une matrice unitaire, le résultat du calcul lui-même ne change pas. $ Q = U_t ^ {-1} ZU_t $.
** 4. Calculez $ QP $ pour calculer le point de contrôle de la nouvelle courbe de Bézier. ** **
\begin{align}
F(t) &= X ~ U_t Q ~P \\
&= X ~ U_t ~ P^{~\prime}
\end{align}
Si $ P ^ {~ \ prime} = QP $, ce sera sous la forme d'une expression représentant la courbe de Bézier. Par conséquent, $ P ^ {~ \ prime} $ est le point de contrôle de la courbe de Bézier après division.
Par la procédure ci-dessus, nous avons pu trouver la ** courbe de Bézier gauche divisée **. Pour trouver la courbe de Bézier sur le côté droit, séparez $ z $ à l'étape 2 comme suit. Trouvez la bonne matrice $ Q ^ {~ \ prime} $.
\begin{align}
F(t) &= \left(
\begin{array}{ccc}
1 & t & t^2 & \dots & t^n
\end{array}
\right)U_tP \\ \\
&= \left(
\begin{array}{ccc}
1 & (z + (1-z)\cdot t) & (z + (1-z)\cdot t)^2 & \dots & (z + (1-z)\cdot t)^n
\end{array}
\right)U_tP
\end{align}
Cependant, vous pouvez réellement trouver $ Q ^ {~ \ prime} $ en utilisant le $ Q $ calculé sans avoir à le calculer directement. En particulier,
Faites simplement l'opération. Comme il est difficile d'expliquer à l'aide de formules mathématiques, reportez-vous à la section de l'exemple de calcul suivant.
Pour ceux qui disent, "$ n $ je ne comprends pas seulement l'histoire suivante!", Voici un exemple de division par calcul matriciel.
Cette fois, à titre d'exemple
Tout d'abord, cette courbe de Bézier est représentée sous la forme d'une matrice, et le point de division $ z = 0,3 $ est séparé.
\begin{align}
F(t) &= XZU_tP \\
&=
\left(
\begin{array}{ccc}
1 & t & t^2
\end{array}
\right)
\left(
\begin{array}{ccc}
1 & 0 & 0 \\
0 & z & 0 \\
0 & 0 & z^2
\end{array}
\right)
\left(
\begin{array}{ccc}
{2 \choose 0}{2 \choose 0}(-1)^0 & 0 & 0 \\
{2 \choose 0}{2 \choose 1}(-1)^1 & {2 \choose 1}{1 \choose 0}(-1)^0 & 0 \\
{2 \choose 0}{2 \choose 2}(-1)^2 & {2 \choose 1}{1 \choose 1}(-1)^1 & {2 \choose 2}{0 \choose 0}(-1)^0
\end{array}
\right)
\left(
\begin{array}{ccc}
P_0 \\
P_1 \\
P_2
\end{array}
\right) \\ \\
&=
\left(
\begin{array}{ccc}
1 & t & t^2
\end{array}
\right)
\left(
\begin{array}{ccc}
1 & 0 & 0 \\
0 & 0.3 & 0 \\
0 & 0 & 0.09
\end{array}
\right)
\left(
\begin{array}{ccc}
1 & 0 & 0 \\
-2 & 2 & 0 \\
1 & -2 & 1
\end{array}
\right)
\left(
\begin{array}{ccc}
P_0 \\
P_1 \\
P_2
\end{array}
\right)
\end{align}
Calculez maintenant $ Q = U_t ^ {-1} ZU_t $.
\begin{align}
Q &= U_t^{-1}ZU_t \\ \\
&=
\left(
\begin{array}{ccc}
1 & 0 & 0 \\
1 & \frac{1}{2} & 0 \\
1 & 1 & 1
\end{array}
\right)
\left(
\begin{array}{ccc}
1 & 0 & 0 \\
0 & 0.3 & 0 \\
0 & 0 & 0.09
\end{array}
\right)
\left(
\begin{array}{ccc}
1 & 0 & 0 \\
-2 & 2 & 0 \\
1 & -2 & 1
\end{array}
\right) \\ \\
&=
\left(
\begin{array}{ccc}
1 & 0 & 0 \\
0.7 & 0.3 & 0 \\
0.49 & 0.42 & 0.09
\end{array}
\right)
\end{align}
Transformez la formule pour en faire une formule en utilisant $ Q $.
\begin{align}
F(t) &= XZU_tP \\
&= XU_tU_t^{-1}ZU_tP \\
&= XU_tQP \\
&=
\left(
\begin{array}{ccc}
1 & t & t^2
\end{array}
\right)
\left(
\begin{array}{ccc}
1 & 0 & 0 \\
-2 & 2 & 0 \\
1 & -2 & 1
\end{array}
\right)
\left(
\begin{array}{ccc}
1 & 0 & 0 \\
0.7 & 0.3 & 0 \\
0.49 & 0.42 & 0.09
\end{array}
\right)
\left(
\begin{array}{ccc}
P_0 \\
P_1 \\
P_2
\end{array}
\right)
\end{align}
Si vous calculez ici pour $ QP $,
\begin{align}
QP &=
\left(
\begin{array}{ccc}
1 & 0 & 0 \\
0.7 & 0.3 & 0 \\
0.49 & 0.42 & 0.09
\end{array}
\right)
\left(
\begin{array}{ccc}
P_0 \\
P_1 \\
P_2
\end{array}
\right) \\ \\
&=
\left(
\begin{array}{ccc}
P_0 \\
0.7P_0 + 0.3P_1 \\
0.49P_0 + 0.42P_1 + 0.09P_2
\end{array}
\right) \\ \\
&= P^{~\prime} \\ \\
\end{align}
\begin{align}
P^{~\prime}_x &=
\left(
\begin{array}{ccc}
0 \\
0.7\cdot 0 + 0.3 \cdot 1 \\
0.49\cdot 0 + 0.42\cdot 1 + 0.09\cdot 2
\end{array}
\right)
=
\left(
\begin{array}{ccc}
0 \\
0.3 \\
0.6
\end{array}
\right)
\end{align}
\begin{align}
P^{~\prime}_y &=
\left(
\begin{array}{ccc}
1 \\
0.7\cdot 1 + 0.3 \cdot 4 \\
0.49\cdot 1 + 0.42\cdot 4 + 0.09\cdot 0
\end{array}
\right)
=
\left(
\begin{array}{ccc}
1 \\
1.9 \\
2.17
\end{array}
\right) \\
\end{align}
Ensuite, nous trouverons la courbe de Bézier sur le côté droit. $ U_t ^ {-1} Z ^ {~ \ prime} U_t = Q ^ {~ \ prime} $ pour le côté droit peut être calculé comme suit en utilisant le $ Q $ calculé.
** 1. Emballez l'élément $ Q $ à droite de 0 $. ** **
\begin{align}
Q &=
\left(
\begin{array}{ccc}
1 & 0 & 0 \\
0.7 & 0.3 & 0 \\
0.49 & 0.42 & 0.09
\end{array}
\right)
\Rightarrow
\left(
\begin{array}{ccc}
0 & 0 & 1 \\
0 & 0.7 & 0.3 \\
0.49 & 0.42 & 0.09
\end{array}
\right)
\end{align}
** 2. Retournez les éléments de la matrice emballés vers la droite à l'envers. ** **
\begin{align}
\left(
\begin{array}{ccc}
0 & 0 & 1 \\
0 & 0.7 & 0.3 \\
0.49 & 0.42 & 0.09
\end{array}
\right)
\Rightarrow
\left(
\begin{array}{ccc}
0.49 & 0.42 & 0.09 \\
0 & 0.7 & 0.3 \\
0 & 0 & 1
\end{array}
\right)
= Q^{~\prime}
\end{align}
C'est facile. Vous pouvez utiliser le $ Q ^ {~ \ prime} $ obtenu pour trouver la ** courbe de Bézier droite divisée **.
\begin{align}
Q^{~\prime}P &=
\left(
\begin{array}{ccc}
0.49 & 0.42 & 0.09 \\
0 & 0.7 & 0.3 \\
0 & 0 & 1
\end{array}
\right)
\left(
\begin{array}{ccc}
P_0 \\
P_1 \\
P_2
\end{array}
\right) \\ \\
&=
\left(
\begin{array}{ccc}
0.49P_0 + 0.42P_1 + 0.09P_2 \\
0.7P_1 + 0.3P_2 \\
P_2
\end{array}
\right) \\ \\
&= P^{~\prime\prime} \\ \\
\end{align}
\begin{align}
P^{~\prime\prime}_x &=
\left(
\begin{array}{ccc}
0.49\cdot 0 + 0.42\cdot 1 + 0.09\cdot 2 \\
0.7\cdot 1 + 0.3 \cdot 2 \\
2
\end{array}
\right)
=
\left(
\begin{array}{ccc}
0.6 \\
1.3 \\
2
\end{array}
\right)
\end{align}
\begin{align}
P^{~\prime\prime}_y &=
\left(
\begin{array}{ccc}
0.49\cdot 1 + 0.42\cdot 4 + 0.09\cdot 0 \\
0.7\cdot 4 + 0.3 \cdot 0 \\
0
\end{array}
\right)
=
\left(
\begin{array}{ccc}
2.17 \\
2.8 \\
0
\end{array}
\right)
\end{align}
Puisque $ P_0 = P_0 ^ {~ \ prime} $, $ P_2 ^ {~ \ prime} = P_0 ^ {~ \ prime \ prime} $, $ P_2 = P_2 ^ {~ \ prime \ prime} $ Il semble que vous puissiez bien le diviser. Si vous dessinez la courbe de Bézier obtenue sur le côté gauche et la courbe de Bézier sur le côté droit, ce sera comme indiqué dans la figure ci-dessous. Vous pouvez bien le diviser.
Un exemple de mise en œuvre d'un programme de division à l'aide d'une matrice est présenté ci-dessous.
import numpy as np
from scipy.special import comb
import matplotlib.pyplot as plt
class Bezier:
def __init__(self, points):
points = [np.array(e).flatten() for e in points]
self._n = len(points)
self._points = points
@property
def dims(self):
return self._n - 1
@property
def points(self):
return self._points
@property
def points4matplot(self):
xs, ys = zip(*self.points)
return xs, ys
def _bernstein(self, n, i, t):
return comb(n, i) * t**i * (1 - t)**(n-i)
def bezier_point(self, t):
p = np.zeros(2)
for i, f in enumerate(self.points):
p += self._bernstein(self.dims, i, t) * f
return np.array(p).flatten()
def split_with_matrix(self, t):
M = self.create_Ut()
iM = np.linalg.inv(M)
Z = np.eye(self.dims + 1)
for i in range(self.dims + 1):
Z[i] = Z[i] * t ** i
## Caluclate Q
Q = iM @ Z @ M
xs, ys = self.points4matplot
X = np.array(xs)
Y = np.array(ys)
left_bezier = Bezier(list(zip(Q @ X, Q @ Y)))
## Make Q'
_Q = np.zeros((self.dims + 1, self.dims + 1))
lst = []
for i in reversed(range(self.dims + 1)):
l = [-1] * i + list(range(self.dims + 1 - i))
lst.append(l)
for i, l in enumerate(lst):
mtx = [Q[i][e] if not e == -1 else 0 for e in l]
_Q[i] = np.array(mtx)
_Q = np.flipud(_Q)
right_bezier = Bezier(list(zip(_Q @ X, _Q @ Y)))
return left_bezier, right_bezier
def create_Ut(self, dims=None):
if dims is None:
dims = self.dims
U = np.zeros([dims + 1, dims + 1])
lmbd = lambda n, i, j: comb(n, j) * comb(n-j, i-j) * (-1)**(i - j)
for i in range(dims + 1):
lst = list(range(i+1)) + [-1]*(dims-i)
elems = [lmbd(dims, i, j) if j >= 0 else 0.0 for j in lst]
mtx = np.array(elems)
U[i] = mtx
return U
def plot(self, ax=None, with_ctrl_pt=True, bcolor="black", ccolor="gray", resolution=100):
if ax is None:
_, ax = plt.subplots()
prev_point = None
for t in np.linspace(0, 1, resolution):
bp = self.bezier_point(t)
if prev_point is None:
prev_point = (bp[0], bp[1])
xs, ys = zip(*(prev_point, (bp[0], bp[1])))
ax.plot(xs, ys, '-', color=bcolor)
prev_point = (bp[0], bp[1])
if with_ctrl_pt:
xs, ys = self.points4matplot
ax.plot(xs, ys, 'x--', lw=2, color=ccolor, ms=10)
return ax
def main():
bezier = Bezier([(0, 1), (1, 4), (2, 0)])
ax = bezier.plot()
left_bezier, right_bezier = bezier.split_with_matrix(0.3)
left_bezier.plot(ax, bcolor = "red")
right_bezier.plot(ax, bcolor = "blue")
plt.grid()
plt.show()
if __name__ == "__main__":
main()
Dans cet article, j'ai présenté la division de la courbe de Bézier d'ordre $ n $. Cela n'a peut-être pas été très organisé, mais pardonnez-moi. J'espère que le contenu de cet article sera utile à quelqu'un.
Recommended Posts