【B01】ベジェ曲線

# 【B01】ベジェ曲線
この記事は[アドカレ](https://qiita.com/advent-calendar/2023/hachinohe_kosen)に参加しています。

# ベジェ曲線とは
　以下の画像のような、直線では表現することができない曲線を表現するために、ベジェ曲線というのがよく使われます。
 ![a_care02.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3005402/04dcb935-5971-36df-df5a-36fb9f2ecb21.jpeg)

　ベジェ曲線は制御点というものを設けることで、複雑な曲線も再現することができるので、視覚的にも分かり易いです。

# 一次のベジェ曲線
　まずは一次のベジェ曲線を考えてみましょう。
![a_care03.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3005402/3b0b4ac6-7b13-0db3-eda0-e8bccf45c682.jpeg)

　一次のベジェ曲線は、`p0`から`p1`までを補間する一本の直線です。
```math
p=f(t)　　(0 \leqq t \leqq 1)
```
引数を、0~1の値をとる`t`として、関数`f(t)`を求めてみます。

`t=0`のとき`f(0)=p0`、`t=1`のとき`f(1)=p1`になると`p0`と`p1`を補間したことになります。これは、一次の関数として考えることができます。一次関数の係数を`a`、`b`として、`f(t)`は
```math
p=f(t)=at+b
```
…のように書くことができます。

このとき、`t=0`で`f(0)=p0`なので、
```math
\begin{align}
p=f(t)&=at+b=p0 \\
&=a \times 0+b=p0 \\
&=b=p0
\end{align}
```
よって、`b=p0`が成り立ちます。

　
次に、`t=1`で`f(1)=p1`なので、
```math
\begin{align}
p=f(t)&=at+b=p1 \\
&=a \times 1+p0=p1 \\
&=a+p0=p1
\end{align}
```
`a+p0=p1`について解くと、
```math
\begin{align}
a+p0&=p1 \\
a&=p1-p0
\end{align}
```
よって、`a=p1-p0`となります。

一次関数の係数である`a`と`b`が求まったので、`p=f(t)`を`t`、`p0`、`p1`で表わすことができます。
```math
\begin{align}
p=f(t)&=at+b \\
&=(p1-p0)t+p0
\end{align}
```

これで、一次のベジェ曲線の式を求めることができました。
```math
p=(p1-p0)t+p0
```
`p0`から`p1`までを`t`(0~1)で補間する関数です。

ここでは二次以降のベジェ曲線で扱いやすくするために、`p=(p1-p0)t+p0`を変形します。
```math
\begin{align}
p&=(p1-p0)t+p0 \\
p&=p1 \times t -p0 \times t + p0 \\
p&=p1 \times t +p0 -p0 \times t \\
p&=p1 \times t +(1-t)p0 \\
p&=t \times p1+(1-t)p0 \\
p&=(1-t)p0+t \times p1
\end{align}
```

# 二次のベジェ曲線
　一次のベジェ曲線の式を使って、二次のベジェ曲線を考えてみましょう。二次のベジェ曲線は、`p0`から`p2`を補完する曲線であり、制御点として`p1`があります。
![a_care04.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3005402/a7920ab6-5038-9bac-57aa-28d952178fad.jpeg)

曲線は、制御点である`p1`を必ず通る訳ではありませんが、`p1`の位置によって曲線の形が変わります。

さて、`p0`と`p2`の補間を考えるには、まず`p0`と`p1`の補間を考えます。
`t=0`のとき`f(t)=p0`、`t=1`のとき`f(t)=p1`となる関数は、一次ベジェ曲線の式を使用して、
```math
p01=(1-t)p0+t \times p1　　・・・(1)
```
と表すことができます。

同様にして、`p1`と`p2`を補間する関数を考えます。
```math
p12=(1-t)p1+t \times p2　　・・・(2)
```

`p0`と`p1`の補間である`p01`と、`p1`と`p2`の補間である`p12`を更に補間します。
```math
p=(1-t)p01+t \times p12　　・・・(3)
```
(1)式、(2)式より、
```math
\begin{align}
p&=(1-t)p01+t \times p12 \\
p&=(1-t)((1-t)p0+t \times p1)+t \times ((1-t)p1+t \times p2) \\
p&=(1-t)^2p0+t(1-t)p1+t(1-t)p1+t^2p2 \\
p&=(1-t)^2p0+2t(1-t)p1+t^2p2
\end{align}
```

これが二次のベジェ曲線の式です。
![a_care05.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3005402/3edba311-42c1-7edd-5c05-75418fafc863.jpeg)
図で見てみると、制御点`p1`の役割がよく分かると思います。

`p0`と`p1`の補間である`p01`は橙色です。
`p1`と`p2`の補間である`p12`は水色です。
`p01`と`p12`の補間である`p`は緑色です。

まだビジュアル的に理解できていない方は、以下のリンク先にあるアニメーションやスライダーを使用して、`t`の値の変化によってベジェ曲線がどのような軌跡を描くのか確認してみてください。
[A Primer on Bézier Curves](https://pomax.github.io/bezierinfo/)
[ベジェ曲線](https://ja.javascript.info/bezier-curve)←このサイトのアニメーションがおすすめです。

# 三次のベジェ曲線
　三次のベジェ曲線は`p0`と`p3`を補間します。制御点は`p1`と`p2`の２つです。
![a_care06.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3005402/51711931-3e32-6645-168c-e9a0160f7ddd.jpeg)

二次のベジェ曲線をビジュアル的に理解していれば、三次はその応用です。

`p0`と`p1`の補間である`p01`は紫色です。
`p1`と`p2`の補間である`p12`は赤色です。
`p2`と`p3`の補間である`p23`は桃色です。
`p01`と`p12`の補間である`p012`は橙色です。
`p12`と`p23`の補間である`p123`は水色です。
`p012`と`p123`の補間である`p`は緑色です。

数式でもみてみましょう。

まず、始点`p0`と制御点`p1`を補間します。
```math
p01=(1-t)p0+t \times p1　　・・・(1)
```

制御点`p1`と制御点`p2`を補間します。
```math
p12=(1-t)p1+t \times p2　　・・・(2)
```

制御点`p2`と終点`p3`を補間します。
```math
p23=(1-t)p2+t \times p3　　・・・(3)
```
　
　
はい。`p01`、`p12`、`p23`を求めることができました。更にこれらを補間します。

`p01`と`p12`を補間します。
```math
p012=(1-t)p01+t \times p12　　・・・(4)
```

`p12`と`p23`を補間します。
```math
p123=(1-t)p12+t \times p23　　・・・(5)
```

最後に、`p012`と`p123`を補完します。
```math
p=(1-t)p012+t \times p123　　・・・(6)
```
(6)式が三次のベジェ曲線の式となります。ですが、このままでは式の全容が分からないので、`t`、`p0`、`p1`、`p2`、`p3`だけで式を表してみましょう。

(4)式、(5)式より、二次のベジェ曲線と同じように変形します。
```math
\begin{align}
p&=(1-t)p012+t \times p123 \\
p&=(1-t)^2p01+2t(1-t)p12+t^2p23
\end{align}
```
(1)式、(2)式、(3)式より、`p01`、`p02`、`p03`にそれぞれ代入していきます。
```math
\begin{align}
p&=(1-t)^2p01+2t(1-t)p12+t^2p23 \\
p&=(1-t)^2((1-t)p0+t \times p1)+2t(1-t)((1-t)p1+t \times p2)+t^2((1-t)p2+t \times p3) \\
p&=(1-t)^3p0+t(1-t)^2p1+2t(1-t)^2p1+2t^2(1-t)p2+t^2(1-t)p2+t^3p3 \\
p&=(1-t)^3p0+3t(1-t)^2p1+3t^2(1-t)p2+t^3p3
\end{align}
```

三次のベジェ曲線の式を求めることができました。

ベジェ曲線は次数を上げていけば上げるほど複雑な曲線を描くことができますが、実際には計算量が少ないながらも充分な自由度を持つ三次の曲線がよく使われます。

# n次のベジェ曲線
　今までは一次から三次のベジェ曲線の式についてみてきました。
```math
\begin{align}
p&=(1-t)p0+t \times p1 \\
p&=(1-t)^2p0+2t(1-t)p1+t^2p2 \\
p&=(1-t)^3p0+3t(1-t)^2p1+3t^2(1-t)p2+t^3p3
\end{align}
```

三つの式を眺めてみると、なにかしらの規則性がありそうなことに気が付くと思います。

n次のベジェ曲線についてシグマを使用して書いてみると、
```math
 \begin{align}
B_{i}&=t^i(1-t)^{n-i} \\
{}_nC_i&=\frac{n!}{i!(n-i)!}
\end{align}
```
として、
```math
p=\sum_{i=0}^{n} {}_nC_i B_{i}p_{i}
```
のようになります。

${}_nC_i$は[二項係数](https://fussy.web.fc2.com/algo/stat1_combination.htm#3)であり、$CB$の部分を[バーンスタイン多項式](https://ja.wikipedia.org/wiki/バーンスタイン多項式)といいます。

# 行列でみるベジェ曲線
　三次のベジェ曲線の式を行列を使って表してみたいと思います。
```math
p=(1-t)^3p0+3t(1-t)^2p1+3t^2(1-t)p2+t^3p3
```
まずは、係数や`t`がごちゃごちゃしていて分かりづらいので、
```math
p=at^3+bt^2+ct+d
```
のような三次関数の式を目指して式を整理します。

`(1-t)`のあたりを展開して、式整理します。 
```math
\begin{align}
p&=(1-t)^3p0+3t(1-t)^2p1+3t^2(1-t)p2+t^3p3 \\
p&=(-t^3+3t^2-3t+1)p0+3t(t^2-2t+1)p1+3t^2(-t+1)p2+t^3p3 \\
p&=(-t^3+3t^2-3t+1)p0+(3t^3-6t^2+3t)p1+(-3t^3+3t^2)p2+(t^3)p3
\end{align}
```
`p0`から`p3`の各係数では、$at^3+bt^2+ct+d$のような三次関数の式になっているので、以下のように書き直せます。
```math
\begin{align}
p&=(-t^3+3t^2-3t+1)p0+(3t^3-6t^2+3t)p1+(-3t^3+3t^2)p2+(t^3)p3 \\
p&=(-p0+3p1-3p2+p3)t^3+(3p0-6p1+3p2)t^2+(-3p0+3p1)t+(p0)
\end{align}
```
 
ここで、以下のような行列の特性を用います。
```math
ap0+bp1+cp2+dp3=

\begin{bmatrix}
a & b & c & d \\
\end{bmatrix}

\begin{bmatrix}
p0 \\ p1 \\ p2 \\ p3
\end{bmatrix}
```
適用してみると、
```math
\begin{align}
p&=(-p0+3p1-3p2+p4)t^3+(3p0-6p1+3p2)t^2+(-3p0+3p1)t+(p0) \\
p&=(\begin{bmatrix} -1 & 3 & -3 & 1 \end{bmatrix}
\begin{bmatrix} p0 \\ p1 \\ p2 \\ p3 \end{bmatrix})t^3+
(\begin{bmatrix} 3 & -6 & 3 & 0 \end{bmatrix}
\begin{bmatrix} p0 \\ p1 \\ p2 \\ p3 \end{bmatrix})t^2+
(\begin{bmatrix} -3 & 3 & 0 & 0 \end{bmatrix}
\begin{bmatrix} p0 \\ p1 \\ p2 \\ p3 \end{bmatrix})t+
(\begin{bmatrix} 1 & 0 & 0 & 0 \end{bmatrix}
\begin{bmatrix} p0 \\ p1 \\ p2 \\ p3 \end{bmatrix})
\end{align}
```
更に適用してみると、
```math
\begin{align}
p&=\begin{bmatrix} t^3 & t^2 & t & 1 \end{bmatrix}
\begin{bmatrix}
\begin{bmatrix} -1 & 3 & -3 & 1 \end{bmatrix}
\begin{bmatrix} p0 \\ p1 \\ p2 \\ p3 \end{bmatrix}
\\
\begin{bmatrix} 3 & -6 & 3 & 0 \end{bmatrix}
\begin{bmatrix} p0 \\ p1 \\ p2 \\ p3 \end{bmatrix}
\\
\begin{bmatrix} -3 & 3 & 0 & 0 \end{bmatrix}
\begin{bmatrix} p0 \\ p1 \\ p2 \\ p3 \end{bmatrix}
\\
\begin{bmatrix} 1 & 0 & 0 & 0 \end{bmatrix}
\begin{bmatrix} p0 \\ p1 \\ p2 \\ p3 \end{bmatrix}
\end{bmatrix} \\

p&=\begin{bmatrix} t^3 & t^2 & t & 1 \end{bmatrix}
\begin{bmatrix}
-1 & 3 & -3 & 1 \\
3 & -6 & 3 & 0 \\
-3 & 3 & 0 & 0 \\
1 & 0 & 0 & 0 
\end{bmatrix}
\begin{bmatrix} p0 \\ p1 \\ p2 \\ p3 \end{bmatrix}

\end{align}
```
綺麗な見た目になりました。ここで以下ように置き換えます。
```math
\begin{align}
T&=\begin{bmatrix} t^3 & t^2 & t & 1 \end{bmatrix} \\
A&=\begin{bmatrix}
-1 & 3 & -3 & 1 \\
3 & -6 & 3 & 0 \\
-3 & 3 & 0 & 0 \\
1 & 0 & 0 & 0 
\end{bmatrix} \\
P&=\begin{bmatrix} p0 \\ p1 \\ p2 \\ p3 \end{bmatrix}
\end{align}
```

最終的に以下のようになります。
```math
p=TAP
```
# 制御点から三次関数の係数
　4つの制御点`p0`~`p3`の値から、三次関数の係数`a`~`d`を求めてみます。

 三次のベジェ曲線の式は以下のように表せるのでした。
```math
\begin{align}
p&=TAP
\end{align}
```
ここで行列を意識してみると、
```math
\begin{align}
at^3+bt^2+ct+d&=TAP \\

\begin{bmatrix} t^3 & t^2 & t & 1 \\ \end{bmatrix}
\begin{bmatrix} a \\ b \\ c \\ d \end{bmatrix} &=
\begin{bmatrix} t^3 & t^2 & t & 1 \end{bmatrix}
\begin{bmatrix}
-1 & 3 & -3 & 1 \\
3 & -6 & 3 & 0 \\
-3 & 3 & 0 & 0 \\
1 & 0 & 0 & 0 
\end{bmatrix}
\begin{bmatrix} p0 \\ p1 \\ p2 \\ p3
\end{bmatrix} \\

\begin{bmatrix} t^3 & t^2 & t & 1 \\ \end{bmatrix}
\begin{bmatrix} a \\ b \\ c \\ d \end{bmatrix} &=
\begin{bmatrix} t^3 & t^2 & t & 1 \end{bmatrix}
\begin{bmatrix}
-p0 & 3p1 & -3p2 & p3 \\
3p0 & -6p1 & 3p2 & 0 \\
-3p0 & 3p1 & 0 & 0 \\
p0 & 0 & 0 & 0 
\end{bmatrix}
\end{align}
```
よって、三次関数の係数は以下のように表すことができます。
```math
\begin{align}
a&=-p0 + 3p1 -3p2 + p3 \\
b&= 3p0 -6p1 + 3p2 \\
c&= -3p0 + 3p1 \\
d&= p0
\end{align}
```

# 三次関数の係数から制御点
　三次関数の係数`a`~`d`から、4つの制御点`p0`~`p3`の値を求めてみます。
```math
\begin{align}
p&=TAP \\
at^3+bt^2+ct+d&=TAP \\

\begin{bmatrix} t^3 & t^2 & t & 1 \\ \end{bmatrix}
\begin{bmatrix} a \\ b \\ c \\ d \end{bmatrix} &=TAP \\
T \begin{bmatrix} a \\ b \\ c \\ d \end{bmatrix} &=TAP \\
\begin{bmatrix} a \\ b \\ c \\ d \end{bmatrix} &=AP
\end{align}
```
 $A$の逆行列$A^{-1}$を両辺に掛けて、
```math
\begin{align}
A^{-1} \begin{bmatrix} a \\ b \\ c \\ d \end{bmatrix} &=A^{-1} AP \\
A^{-1} \begin{bmatrix} a \\ b \\ c \\ d \end{bmatrix} &=P \\
\end{align}
```
よって、三次関数の係数`a`~`d`から、4つの制御点`p0`~`p3`の値を求めるには以下のような式を用いります。
```math
\begin{align}
P&=A^{-1} \begin{bmatrix} a \\ b \\ c \\ d \end{bmatrix}\\
\begin{bmatrix} p0 \\ p1 \\ p2 \\ p3 \end{bmatrix} &=
\begin{bmatrix}
0 & 0 & 0 & 1 \\
0 & 0 & \frac{1}{3} & 1 \\
0 & \frac{1}{3} & \frac{2}{3} & 1 \\
1 & 1 & 1 & 1 
\end{bmatrix}
\begin{bmatrix} a \\ b \\ c \\ d \end{bmatrix}
\end{align}
```

# ベジェを座標系に落とし込む
　n次のベジェ曲線の式は以下のようになるのでした。
```math
p=\sum_{i=0}^{n} {}_nC_i B_{i}p_{i}
```
ですが、これは一次元の座標系での式となります。

二次元平面でのベジェ曲線は以下の式で表すことができます。
```math
\begin{bmatrix} x \\ y \end{bmatrix}=
\begin{bmatrix}
\sum_{i=0}^{n} {}_nC_i B_{i}x_{i} \\
\sum_{i=0}^{n} {}_nC_i B_{i}y_{i} \end{bmatrix}
```

同様にして、三次元空間では以下のようになります。
```math
\begin{bmatrix} x \\ y \\ z \end{bmatrix}=
\begin{bmatrix}
\sum_{i=0}^{n} {}_nC_i B_{i}x_{i} \\
\sum_{i=0}^{n} {}_nC_i B_{i}y_{i} \\
\sum_{i=0}^{n} {}_nC_i B_{i}z_{i} \end{bmatrix}
```

# プログラム例
　三次のベジェ曲線関連のプログラム例を示します。

まずは制御点$P$を用いるベジェの関数です。もとの式のままだと掛け算を14回、足し算を3回するので、それよりも演算の回数が少なくなるように式変形します。
```math
\begin{align}
p&=(1-t)^3p0+3t(1-t)^2p1+3t^2(1-t)p2+t^3p3 \\
p&=(1-t)^3p0+t(3(1-t)^2p1+3t(1-t)p2+t^2p3) \\
p&=(1-t)^3p0+t(3(1-t)^2p1+t(3(1-t)p2+tp3))
\end{align}
```
足し算は3回のままですが、掛け算を11回まで減らすことができました。

プログラムでは、引数`t`の値に対応するベジェ曲線上の位置を返します。
```c++
//return p[0]*s^3 + t(3*p[1]*s^2 + t(3*p[2]*s + t*p[3]))
double BezierLineP(double* p, double t) {
	double s = 1 - t;
	return *(p + 0) * s * s * s + t * (3 * *(p + 1) * s * s + t * (3 * *(p + 2) * s + t * *(p + 3)));
}
```

次は係数`a`~`d`を用いた三次関数です。もとの式のままだと掛け算を6回、足し算を3回するので、それよりも演算の回数が少なくなるように式変形します。
```math
\begin{align}
p&=at^3+bt^2+ct+d \\
p&=(at^2+bt+c)t+d \\
p&=((at+b)t+c)t+d \\
\end{align}
```
足し算は3回のままですが、掛け算を3回まで減らすことができました。

```c++
//return ((at+b)t+c)t+d
double BezierLineC(double* c, double t) {
	return ((*(c + 0) * t + *(c + 1)) * t + *(c + 2)) * t + *(c + 3);
}
```

最後は、制御点$P$から三次関数の係数を求める関数です。
```c++
//p->c
void BezierCoefficient(double* p, double* c) {
	*(c + 0) = -1 * *(p + 0) + 3 * *(p + 1) - 3 * *(p + 2) + 1 * *(p + 3);
	*(c + 1) = 3 * *(p + 0) - 6 * *(p + 1) + 3 * *(p + 2);
	*(c + 2) = -3 * *(p + 0) + 3 * *(p + 1);
	*(c + 3) = 1 * *(p + 0);
}
```
制御点$P$からベジェ曲線を求めることもできますが、計算回数は三次関数の係数を用いた方法のほうが少ないです。場合によって使い分けるのが大事です。

# 参考文献
[ベジェ曲線 ウィキペディア](https://ja.wikipedia.org/wiki/ベジェ曲線)
[曲線描画(2) ベジェ曲線](https://fussy.web.fc2.com/algo/curve2_bezier.htm)
[A Primer on Bézier Curves](https://pomax.github.io/bezierinfo/)
[一から学ぶベジェ曲線](https://postd.cc/bezier-curves/)
[ベジェ曲線](https://ja.javascript.info/bezier-curve)

# むすび
　ベジェ曲線について解説しました。今更ですが、図はイメージです。実際のベジェの世界はもっとすごいので、興味のある方は参考文献を読んでみてください。