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. 図0.png

** 1. Find the point that is $ t: 1-t $ on the line connecting the control points of the Bezier curve (point $ Q $). ** ** 図1.png

** 2. For the line connecting the obtained points, find the point that is $ t: 1-t $ (point $ R $). ** ** 図2.png

** 3. Repeat until there is one point (point $ S $). ** ** 図3.png

** 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. ** ** 図4.png

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. 図5.png

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. 図6.png

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,

  1. Pack the elements of $ Q $ to the right by $ 0 $.
  2. 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 P_0=(0, 1),~ P_1=(1, 4),~ P_2=(2, 0) Suppose you want to divide a quadratic Bezier curve with a control point at $ t = 0.3 $. It is a curve as shown in the figure below.

f1.png

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.

f2.png

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

Recommended Posts

I want to divide an nth-order Bezier curve at an arbitrary point ~ Recursion and matrix ~
I want to find the intersection of a Bezier curve and a straight line (Bezier Clipping method)
I want to leave an arbitrary command in the command history of Shell
I want to make an automation program!
I want to make a music player and file music at the same time
I want to write an element to a file with numpy and check it.
I want to have recursion come to my mind
I wanted to play with the Bezier curve
I want to be an OREMO with setParam!
I want to convert an image to WebP with lollipop
I want to handle optimization with python and cplex
I want to develop an Android application on Android (debugging)
I want to publish the product at the lowest cost
I want to improve efficiency with Python even in an experimental system (2) RS232C and pySerial
I want to copy an English paper from pdf and put it in Google Translate