3次ベジェ曲線の接続 - 機械学習基礎理論独習

要件

制御点 ${\bf p}_0,{\bf p}_1,{\bf p}_2,{\bf p}_3\in{\mathbb R}^2$ からなる3次ベジェ曲線の名前を $\rm P$ 、位置ベクトルを ${\bf p}(t)$ 、
制御点 ${\bf q}_0,{\bf q}_1,{\bf q}_2,{\bf q}_3\in{\mathbb R}^2$ からなる3次ベジェ曲線の名前を $\rm Q$ 、位置ベクトルを ${\bf q}(u)$ とします。
${\bf p}_0,{\bf p}_1,{\bf p}_2,{\bf p}_3$ を更新し ${\bf p}_3={\bf q}_0$ となるように曲線を接続することを考えます。

接続の種類

C0連続:位置ベクトルの0階微分連続
C1連続:位置ベクトルの0,1階微分が連続
C2連続:位置ベクトルの0,1,2階微分が連続
G0連続:C0連続と同じ
G1連続:接ベクトルが同じ方向
G2連続:符号付曲率が連続

G1連続はG0連続であり、G2連続はG0,G1連続でもあるとします。
C0連続とG0連続は同じことだと思うので、以下G0連続という言葉は使わず、C0連続という言葉を使います。

C0連続

$\rm P$ の終点と $\rm Q$ の始点が一致するように $\rm Q$ 全体を動かします。
※ $\rm Q$ 全体を動かす必要は無いのですが、全体を動かした方が使いやすいと思うので。
${\bf p}(1)={\bf p}_3,{\bf q}(0)={\bf q}_0$ なので、 $\Delta={\bf p}_3-{\bf q}_0$ とし、 ${\bf q}_i={\bf q}_i+\Delta\ (i=0,1,2,3)$ とすればよいです。

C1連続

まず、C0連続にします。
次に ${\bf q}_1$ を動かして1階微分が一致するようにします。
$\dfrac{d{\bf p}}{dt}(1)=\dfrac{d{\bf q}}{du}(0)$ を ${\bf q}_1$ について解きます。

$\begin{eqnarray} &&\dfrac{d{\bf p}}{dt}(1)=\dfrac{d{\bf q}}{du}(0)\\ &&\Rightarrow 3({\bf p}_3-{\bf p}_2)=3({\bf q}_1-{\bf q}_0)\\ &&\Rightarrow{\bf q}_1=2{\bf p}_3-{\bf p}_2\tag{1} \end{eqnarray}$

式 $(1)$ のように ${\bf q}_1$ を更新すればよいです。

C2連続

まず、C0,C1連続にします。
次に ${\bf q}_2$ を動かして2階微分が一致するようにします。
$\dfrac{d^2{\bf p}}{dt^2}(1)=\dfrac{d^2{\bf q}}{du^2}(0)$ を ${\bf q}_2$ について解きます。

$\begin{eqnarray} &&\dfrac{d^2{\bf p}}{dt^2}(1)=\dfrac{d^2{\bf q}}{du^2}(0)\\ &&\Rightarrow 6 ( ({\bf p}_3-{\bf p}_2)-({\bf p}_2-{\bf p}_1 ) ) = 6 ( ({\bf q}_2-{\bf q}_1)-({\bf q}_1-{\bf q}_0 ) )\\ &&\Rightarrow{\bf q}_2=4{\bf p}_3-4{\bf p}_2 + {\bf p}_1\tag{2} \end{eqnarray}$

式 $(2)$ のように ${\bf q}_2$ を更新すればよいです。

G1連続

まず、C0連続にします。
次に ${\bf q}_1$ を動かして接ベクトルの方向が一致するようにします。
${\bf q}_0$ を中心に ${\bf q}_1$ を回転させて、 ${\bf q}_1-{\bf q}_0$ が ${\bf p}_3-{\bf p}_2$ と同じ方向になるようにします。
${\bf p}_{23}={\bf p}_3-{\bf p}_2,{\bf q}_{01}={\bf q}_1-{\bf q}_0$ とおくと、
${\rm atan2}({\bf p}_{23y},{\bf p}_{23x})-{\rm atan2}({\bf q}_{01y},{\bf q}_{01x})$ が回転角度です。

G2連続

※G2連続では3つのアルゴリズムを紹介します。
　「アルゴリズムQ1」は条件付きで ${\bf q}_1$ のみを動かすアルゴリズムです。
　「アルゴリズムQ2」は ${\bf q}_2$ のみを動かすアルゴリズムです。
　「アルゴリズムQ1Q2」は ${\bf q}_1,{\bf q}_2$ の両方を動かすアルゴリズムです。

まず、C0,G1連続にします。
次に ${\bf q}_1, {\bf q}_2$ を動かして曲率が一致するようにします。
${\bf q}_0$ を原点、 ${\bf q}_1$ が x+ 上に来るように ${\bf q}_i\ (i=0,1,2,3)$ を座標変換します。
座標変換の平行移動は $-{\bf q}_0$ であり、回転角度 $\theta$ は ${\bf q}_{01}={\bf q}_1-{\bf q}_0$ とおくと、 $\theta=-{\rm atan2}({\bf q}_{01y},{\bf q}_{01x})$ です。
よって、座標変換を表す行列 ${\bf M}$ は

$\begin{eqnarray} {\bf M}=\begin{pmatrix} \cos\theta & -\sin\theta & 0\\ \sin\theta & \cos\theta & 0\\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} 0 & 0 & -p_{0x}\\ 0 & 0 & -p_{0y}\\ 0 & 0 & 1 \end{pmatrix}\tag{3} \end{eqnarray}$

となります。
下図は座標変換後の曲線 $\rm Q$ です。
座標変換後の ${\bf q}_1,{\bf q}_2$ の座標を ${\bf q}_1=(x, 0),{\bf q}_2=(q_{2x}, y)$ とおきます。
${\bf p}_3$ における符号付曲率半径を $R$ とします。
$R=\dfrac{\left(x'^2+y'^2\right)^{\frac{3}{2}}}{x'y''-y'x''}$ であり、変換後の座標系では、 $x'=3x,y'=0,y''=6y$ なので、
$R=\dfrac{3x^2}{2y}$ となります。

$y$ と $R$ が同じ符号の時、 ${\bf q}_1$ を更新します。(アルゴリズムQ1と命名します)
$\hat{x}=\sqrt{\dfrac{2}{3}yR}$ となるので、求める $\hat{\bf q}_1$ は変換後の座標系で $\hat{\bf q}_1=(\sqrt{\dfrac{2}{3}yR}, 0)$ となります。
これをもとの座標系にすると、 ${\bf M}^{-1}\hat{\bf q}_1$ となります。

$y$ と $R$ が異なる符号の時、 ${\bf q}_2$ を更新します。(アルゴリズムQ2と命名します)
※アルゴリズムQ2は $y$ と $R$ の符号は関係ないアルゴリズムです。
　あなたが ${\bf q}_1$ を更新したいのであれば、最初からアルゴリズムQ2を採用しましょう。
　 ${\bf q}_1$ を動かせないときに、やむを得ず ${\bf q}_2$ を動かすという風に紹介したかったので、「異なる符号」と書いています。
$\hat{y}=\dfrac{3x^2}{2R}$ となるので、求める $\hat{\bf q}_2$ は変換後の座標系で $\hat{\bf q}_2=(q_{2x}, \dfrac{3x^2}{2R})$ となります。
これをもとの座標系にすると、 ${\bf M}^{-1}\hat{\bf q}_2$ となります。

$y$ と $R$ の関係性により、 ${\bf q}_1$ を更新したり、 ${\bf q}_2$ を更新するというアルゴリズムですが
これだと場合によっては、 ${\bf q}_1$ または ${\bf q}_2$ が大きく更新されて曲線の形状が大きく変化する場合があります。
それを避けるために ${\bf q}_1,{\bf q}_2$ の片方ではなく、両方更新することによって、
曲線の形状をあまり変化させないようにG2連続してみたいと思います。(アルゴリズムQ1Q2と命名します)
$\hat{\bf q}_1=(x+\Delta x, 0), \hat{\bf q}_2=(q_{2x},y+\Delta y)$ とします。
$x+\Delta x\geq 0$ より、 $\Delta x\geq x$ であることに注意しましょう。 $\Delta y$ は拘束条件はありません。

目的関数 ${\rm E}$ を以下のように定義します。

$\begin{eqnarray} {\rm E}=\Delta x^2+\Delta y^2\tag{4} \end{eqnarray}$

また、曲率の式より、 $\Delta y$ を $\Delta x$ で表せます。

$\begin{eqnarray} &&R=\dfrac{3(x+\Delta x)^2}{2(y+\Delta y)}\\ &&\Rightarrow \Delta y=-y+\dfrac{3(x+\Delta x)^2}{2R}\tag{5} \end{eqnarray}$

式 $(5)$ を式 $(4)$ に代入します。

$\begin{eqnarray} {\rm E}=\frac{9}{4{\rm R}^2}\Delta x^4 + \frac{9x}{{\rm R}^2}\Delta x^3+\left(1-\frac{3y}{R}+\frac{27x^2}{2{\rm R}^2}\right)\Delta x^2+\left(-\frac{6xy}{{\rm R}}+\frac{9x^3}{{\rm R}^2}\right)\Delta x+y^2-\frac{3x^2y}{{\rm R}}+\frac{9x^4}{4{\rm R}^2}\tag{6} \end{eqnarray}$

式 $(6)$ を $\Delta x$ で偏微分して $=0$ とおきます。

$\begin{eqnarray} &&\frac{\partial{\rm E}}{\partial\Delta x}=0\\ &&\Rightarrow\frac{9}{{\rm R}^2}\Delta x^3 + \frac{27x}{{\rm R}^2}\Delta x^2+\left(2-\frac{6y}{\rm R}+\frac{27x^2}{{\rm R}^2}\right)\Delta x-\frac{6xy}{{\rm R}}+\frac{9x^3}{{\rm R}^2}=0\tag{7} \end{eqnarray}$

式 $(7)$ を満たす実数解の内、式 $(6)$ を最小にするものを採用します。
(式 $(7)$ を解くには、解析的解法と数値的解法(ニュートン法など)がありますが、数値的解法の方がいいんじゃないかって思っています。)
式 $(7)$ を解いて $\Delta x$ が定まると、式 $(4)$ より $\Delta y$ も定まります。
そうすると、 $\hat{\bf q}_1,\hat{\bf q}_2$ が定まります。
これらをもとの座標系にすると、 ${\bf M}^{-1}\hat{\bf q}_1,{\bf M}^{-1}\hat{\bf q}_2$ となります。

最後に

・G2連続のアルゴリズムですが個人的には、3次方程式がうまく解けない場合に、
　 ${\bf q}_1,{\bf q}_2$ の片方のみを更新するのが良いのでは思っております。

・G2連続の発展的な話題として、 ${\bf p}_2$ を動かしてみるのもありかもしれません。
　 ${\bf q}_1-{\bf p}_2$ の角度と $|{\bf p}_2-{\bf p}_3|$ と $|{\bf q}_1-{\bf q}_0|$ の3つ自由度による最適化も面白いかもしれない。

・G2連続の発展的な話題として、 ${\bf p}_3={\bf q}_0$ を中心に ${\bf p}_2,{\bf q}_1$ を同じ角度回転させて
　 $y$ と $R$ が同じ符号になるようにして「アルゴリズムQ1」が適用するのはどうかと思いましたが、
　これだと接続部周辺の形状が変化しすぎる場合があり得るので、これはよくないなと思いました。(考えたことをメモしておきます。)

・パソコンのモニタなどで見る曲線の接続はG1連続ぐらいで十分だと思います。

参考リンク

Smoothness - Wikipedia
ベジエ曲線のなめらかさについて社内勉強会で語った
 JavaScriptでニュートン法を用いて三次方程式を解いて、ついでにグラフも描く