[PYTHON] I want to divide an nth-order Bezier curve at an arbitrary point ~ Recursion and matrix ~
Introduction
From certain circumstances, I was investigating the division of the Bezier curve. I couldn't find many sites that explain how to divide a ** $ n $ Bezier curve **.
Well, since Bezier curves are used only up to the 4th order, there is no problem even if you do not know the $ n $ order. I really wanted to implement it neatly, so I investigated it, so I will leave it as knowledge. I have posted an implementation example in the article, but I have also posted it on GitHub (see line_module.py).
Method using recursion
This is a method called De Casteljau's algorithm. This method uses recursion to find the control points for a new Bezier curve after division. This method is ** simple in concept and easy to implement **, but has the disadvantage of being difficult to use in embedded systems with a small stack size **. Well, the Bezier curve is only up to the 4th order (omitted)
algorithm
As an example, let's divide the Bezier curve shown below.
Let $ t $ be the parameter of the Bezier curve at the dividing point, and $ P_i $ be the initial control point of the Bezier curve.
** 1. Find the point that is $ t: 1-t $ on the line connecting the control points of the Bezier curve (point $ Q $). ** **
** 2. For the line connecting the obtained points, find the point that is $ t: 1-t $ (point $ R $). ** **
** 3. Repeat until there is one point (point $ S $). ** **
** 4. At all points including the initial control points, the set of points with the smallest subscript and the set of points with the largest subscript are the control points after division. ** **
For the last part (minimum subscript), it will be easier to understand if you write the following tree structure.
The tree structure shown below shows from which line segment (which two points) a new point was generated.
The outermost points of this tree, namely points $ \ {P_0, Q_0, R_0, S_0 \} $ and points $ \ {S_0, R_1, Q_2, P_3 \} $, become new Bezier curves. I will. You can see that the start and end points of the original Bezier curve are common, and that the point (division point $ S_0 $) in the parameter $ t $ is the start point (end point) of the new Bezier curve.
program
An implementation example of a split program using recursion is shown below.
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()
Method using matrix
This is a method to find a new control point using a matrix. This method ** does not matter if the stack size is small **, but ** the calculation is a little complicated **. Well, the Bezier curve is at most 4th order (omitted)
Represent the Bezier curve in the form of a matrix
As a preliminary step to explain the algorithm, we will explain the matrix representation of Bezier curves. The Bezier curve can be expressed by the following equation.
\begin{align}
F(t) &= \sum_{i=0}^{n} P_i~B_{i}^{n}(t)
\end{align}
$ P_i $ represents the control point and $ B_ {i} ^ {n} (t) $ represents the Bernstein basis function. $ n $ represents the order of the Bezier curve plus $ 1 $ (the number of $ = $ control points). This is expressed in the form of matrix multiplication as follows.
\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}
Also, the Bezier curve formula can be expressed as follows by rearranging it for $ t $. Here, $ a_n $ represents an appropriate coefficient.
\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}
From the above, the following equation can be derived.
\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}
It's pretty clean. Here, we know that $ A $ can be expressed as follows by a lower triangular matrix $ U_t $ and a control point $ 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}
Therefore, finally, the Bezier curve equation can be expressed using a matrix as follows.
F(t) = BP = XU_tP
algorithm
The algorithm for Bezier curve division using a matrix is as follows.
** 1. Arrange the Bezier curve for $ t $ and express it in the form of a matrix. ** **
\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. Separate the division point $ z ~ (z \ int) $. ** **
\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. Transform the formula as follows to find $ 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}
In the above calculation, the position of $ U_t $ is moved. Since $ U_t U_t ^ {-1} $ is an identity matrix, the calculation result itself does not change. $ Q = U_t ^ {-1} ZU_t $.
** 4. Calculate $ QP $ to calculate the control points of the new Bezier curve. ** **
\begin{align}
F(t) &= X ~ U_t Q ~P \\
&= X ~ U_t ~ P^{~\prime}
\end{align}
If $ P ^ {~ \ prime} = QP $, it will be in the form of an expression representing a Bezier curve. Therefore, $ P ^ {~ \ prime} $ is the control point of the Bezier curve after division.
By the above procedure, we were able to find the ** divided left Bezier curve **. To find the Bezier curve on the right side, separate $ z $ in step 2 as follows. Find the matrix $ Q ^ {~ \ prime} $ for the right side.
\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}
However, you can actually find $ Q ^ {~ \ prime} $ using the calculated $ Q $ without having to calculate it directly. In particular,
- Pack the elements of $ Q $ to the right by $ 0 $.
- Flip the elements of the matrix packed to the right upside down.
Just do the operation. Since it is difficult to explain using mathematical formulas, refer to the section of the following calculation example.
Calculation example
For those who say, "$ n $ I don't understand just the next story!", Here is an example of division by matrix calculation.
This time, as an example
First, this Bezier curve is represented in the form of a matrix, and the division point $ z = 0.3 $ is separated.
\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}
Now calculate $ 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}
Transform the formula to make it a formula using $ 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}
If you calculate for $ QP $ here,
\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}
And ** the control point of the Bezier curve on the left side of the division ** $P_0^{~\prime}=(0, 1),~ P_1^{~\prime}=(0.3, 1.9),~ P_2^{~\prime}=(0.6, 2.17)$ It turned out to be.
Next, we will find the Bezier curve on the right side. $ U_t ^ {-1} Z ^ {~ \ prime} U_t = Q ^ {~ \ prime} $ for the right side can be calculated using the calculated $ Q $ as follows.
** 1. Pack the $ Q $ element to the right by $ 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. Flip the elements of the matrix packed to the right upside down. ** **
\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}
It's easy. You can use the calculated $ Q ^ {~ \ prime} $ to obtain the ** divided right Bezier curve **.
\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}
Now the control point of the Bezier curve on the right side is $P_0^{~\prime\prime}=(0.6, 2.17),~ P_1^{~\prime\prime}=(1.3, 2.8),~ P_2^{~\prime\prime}=(2, 0)$ It turned out to be.
Since $ P_0 = P_0 ^ {~ \ prime} $, $ P_2 ^ {~ \ prime} = P_0 ^ {~ \ prime \ prime} $, $ P_2 = P_2 ^ {~ \ prime \ prime} $ It seems that you can divide it well. If you draw the obtained Bezier curve on the left side and the Bezier curve on the right side, it will look like the figure below. You can divide it well.
program
An implementation example of a division program using a matrix is shown below.
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()
in conclusion
In this article, we introduced the division of the $ n $ Bezier curve. It may not have been very organized, but please forgive me. I hope the content of this article is useful to someone.
References
- A Primer on Bezier Curves
- A Matrix method for efficient computation of Bernstein coefficients
- De Casteljau's algorithm (wiki)
Recommended Posts