ベジェ曲線と直線の交点を求めたい（BezierClipping 法）

## はじめに
後輩の研究でベジェ曲線と直線の交点を求める必要がありました．
研究のお手伝いをしていた私は，

「**Pythonなら強々ライブラリとか簡潔な実装例とかあるでしょ！余裕余裕** :grin: 」

とか思っていたのですが，**私にも理解できるレベルで**書いてある記事が見つけられませんでした．
仕方なく，比較的簡単そう＆多少の解説PDFが見つけられたBezierClipping法という手法に焦点を定め，
論文を読み漁って理解・実装したので，その知識を残しておこうと思います．
記事中に多少正確性に欠ける表現があるかもしれませんがご了承ください（そして訂正例を教えて下さい）．

実装したコードは[GitHub](https://github.com/takkaO/Bezier_Clipping_Algorithm)に公開しています．

## ベジェ曲線ってなんぞ？
ベジェ曲線（Bezier Curve）は曲線の表し方の一種です．
アウトラインフォントやCG等コンピュータ上で曲線を表すときによく利用されています．
ベジェ曲線は曲線の始点と終点に加えて，制御点（Control Point）と呼ばれるいくつかの点で構成されます．
以降では始点と終点も含めて制御点と呼称します．
ベジェ曲線は媒介変数 $t$ $(0\leq t \leq 1)$を用いて次のように表されます．

```math
\begin{align}
P(t) &= \sum_{i=0}^{n} F_i B_{i}^{n}(t) \\
&= \sum_{i=0}^{n} F_i {n \choose i} t^i \left(1-t \right)^{n-i}
\end{align}
```

$F_i$ は制御点，$B_{i}^{n}$ はバーンスタイン基底関数（Bernstein Basis Polynomials）です．
また，${n \choose i}$ は組み合わせ（combination）を表します．
世界標準はこの書き方らしいですが，日本の高校で習ったように書くと下のようになります．

```math
P(t) = \sum_{i=0}^{n} F_i ~~ {}_n \mathrm{C} {}_r ~~ t^i \left(1-t \right)^{n-i}
```

ベジェ曲線は $t$ $(0\leq t \leq 1)$における点の集合なので，
$t$ の刻み幅を細かくすればするほど曲線がなめらかになります．
以下に刻み幅を変えた場合の例を出します．灰色の$\times$点は制御点を表しています．
**$(0\leq t \leq 1)$を5分割した場合 (5点だけ計算した場合)．**
![1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/580996/e5166998-d571-49e9-5f50-3cb1c50c85e5.png)

**$(0\leq t \leq 1)$を10分割した場合 (10点だけ計算した場合)．**
![2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/580996/d224b784-251d-03a5-b7b3-f5b3efef381d.png)

**$(0\leq t \leq 1)$を50分割した場合 (50点だけ計算した場合)．**
![3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/580996/cb7bc378-6104-c122-5caf-56ba9937f3e8.png)

**$(0\leq t \leq 1)$を100分割した場合 (100点だけ計算した場合)．**
![4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/580996/5c3f2a83-5a8f-68a5-791d-a5e162c65c73.png)


## BezierClipping 法
### 概要
BezierClipping 法はベジェ曲線の性質を利用して曲線と直線の交点を求める手法です．
西田らによって1990年に提案されました．
特徴としては以下のような点が挙げられています（一部抜粋）．

- 指定した範囲の全ての解（交点）が求出可能
- 安定して解を求められる
- 解の有無の判定が高速
- 解を小さい順に算出可能
- 1次式のみの反復計算　...etc

ベジェ曲線同士の交点や3次元空間における交点などにも応用できるようですが，
今回はベジェ曲線と直線に限定します．

### アルゴリズム
アルゴリズムの大まかな流れは次の通りです．
**1. 対象のベジェ曲線を直線との距離を表すベジェ曲線に変換する（$xy$ 平面から $td$ 平面へ変換する）．**
![5-6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/580996/e7971890-4e55-5eb1-56b6-9608d469db1e.png)
**2. 変換したベジェ曲線の制御点の凸包（convex-hull）を求め，$t$ 軸との交点 $t_{min}, t_{max}$ を求める．**
![7.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/580996/b1053482-8dbb-adeb-c106-bf56b620d23d.png)

**3. 求めた $t_{min}, t_{max}$ でベジェ曲線を分割し，新しいベジェ曲線とする．**
![8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/580996/7d6caefd-a4da-b15a-294a-c426fc99577f.png)
**4. 2-3の処理を区間 $[t_{min},~ t_{max}]$ が十分小さくなるまで繰り返す．**
![9.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/580996/d33c13ba-c434-e91b-9f30-6a4ccb2b470b.png)
**5. 求まった $t$ (上図の赤点) が直線との交点における元のベジェ曲線の媒介変数 $t$ となる．**
![10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/580996/437b1869-74ed-832f-406e-3e2afa875aec.png)

アルゴリズム自体は非常に単純ですね．
おそらく詰まるところは

- **ベジェ曲線の変換処理**
- **複数交点が存在するときの処理**

だと思うので，次節で詳しく説明します．

### ベジェ曲線の変換処理
二次元平面座標空間を前提とします．
平面空間内にある直線 $L$ は次のように表すことができます．$a, b, c$ は定数です．

```math
ax + by + c = 0
```

座標空間内の任意の点 $(x_1, y_1)$ と直線 $L$ の距離 $d$ は次のように表せます．

```math
d = \frac{ax_1+by_1+c}{\sqrt{a^2 + b^2}}
```

ここで，ベジェ曲線上の点 $(x, y)$ は $t$ を媒介変数として次のように求められます．
$(x_i, y_i)$はベジェ曲線の制御点を表します．

```math
\begin{align}
x(t) &= \sum_{i=0}^{n} x_i B_i^n (t) \\
y(t) &= \sum_{i=0}^{n} y_i B_i^n (t)
\end{align}
```

これらを使い，直線 $L$ とベジェ曲線間の距離 $D$ は次のように表せます．

```math
\begin{align}
D &= \frac{ax(t)+by(t)+c}{\sqrt{a^2+b^2}} \\ \\
  &= \frac{1}{\sqrt{a^2+b^2}}\left( a\sum_{i=0}^{n} x_i B_i^n (t) + b\sum_{i=0}^{n} y_i B_i^n (t) + c \right)
\end{align}
``` 

$\sum_{i=0}^{n} B_i^n (t) = 1$ であるので，

```math
\begin{align}
D &= \frac{1}{\sqrt{a^2+b^2}}\left( a\sum_{i=0}^{n} x_i B_i^n (t) + b\sum_{i=0}^{n} y_i B_i^n (t) + c\sum_{i=0}^{n} B_i^n (t) \right) \\ \\
  &= \frac{1}{\sqrt{a^2+b^2}}\sum_{i=0}^{n} (ax_i+by_i+c) B_i^n (t) \\ \\
  &= \frac{1}{\sqrt{a^2+b^2}}\sum_{i=0}^{n} d_i B_i^n (t) 
\end{align}
```

ここで $D$ は制御点 $\left(\frac{i}{n}, d_i \right)$ から成る有理ベジェ曲線となります．
また $D$ は $td$ 平面のベジェ曲線，すなわち横軸 $t$ 縦軸 $d$をとる平面空間におけるベジェ曲線になります．
変換したベジェ曲線 $D$ が，ある $t ~(0 \leq t \leq 1) $ において $d=0$ をとるということは，
その $t$ 値において変換前の $xy$ 平面におけるベジェ曲線と直線間の距離がゼロになることを意味しています．
したがって，直線とベジェ曲線の交点における媒介変数 $t$ が算出できたことになります．
あとは下式に交点における媒介変数 $t$ を入れれば，$xy$ 平面における交点の座標が求まります．

```math
\begin{align}
x(t) &= \sum_{i=0}^{n} x_i B_i^n (t) \\
y(t) &= \sum_{i=0}^{n} y_i B_i^n (t)
\end{align}
```

### 複数点が存在するときの処理
複数の交点が存在するときの処理は簡単です．
複数の交点が存在する場合は $t_{min}$ と $t_{max}$ の変化が小さくなります．
例を下図に示します．図中の青点が1つ前の $t_{min}, t_{max}$ の点を表しています．

**1回目**
![1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/580996/c719cad3-53b8-1fdb-bc9b-9e8e541c1748.png)
**2回目**
![2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/580996/6a434425-d3ef-9129-e652-867e8708e408.png)
**3回目**
![3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/580996/3f355d15-b247-5fd8-3fc6-a00065922780.png)
**4回目**
![4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/580996/5e22aaa9-da69-b0e5-a022-a4ff4e5773b7.png)

3回目から $t_{min}, t_{max}$ の位置がほぼ変わっていないことが分かります．
この場合は曲線を適当な点で分割し，それぞれのベジェ曲線に対して再度 BezierClipping 法を適用します．
分割する点に関しては言及されていないようでした．
私の実装では中点で分割するようにしています．

## おわりに
本記事ではBezierClipping 法について解説しました．
私は最初の変換の部分がなかなか理解できずにかなり苦労しました．
本記事の内容が誰かのお役に立てれば幸いです．


## 参考文献
- [Bezier Clipping](http://nishitalab.org/user/nis/ourworks/BezClip/BezierClipping.html)
- [Bezier Clipping およびその計算機支援形状設計への応用](https://www.ieice.org/jpn/event/FIT/2016/data/pdf/I-011.pdf)
- [Curve intersection using Bezier clipping](http://nishitalab.org/user/nis/cdrom/cad/CAGD90Curve.pdf)
- [曲線・曲面の処理](https://www.jstage.jst.go.jp/article/bjsiam/13/3/13_KJ00003509912/_pdf)
- [BEZIER CURVES INTERSECTION USING MONGE PROJECTION](http://www.sccg.sk/~kg/palaj/pub/3_Kocovce%202007.pdf)
- [ベジェ曲線を手で描いてみる - にせねこメモ](https://nixeneko.hatenablog.com/entry/2015/06/26/075022)
- [A Primer on Bezier Curves](https://pomax.github.io/bezierinfo/#splitting)
- [A Matrix method for efficient computation of Bernstein coefficients](https://interval.louisiana.edu/reliable-computing-journal/volume-17/reliable-computing-17-pp-40-71.pdf)