点列に3次ベジェ曲線をフィットさせる

はじめに

本アルゴリズムは、何か文献を参考にしたものではないので、間違っている可能性があります。予めご了承ください。

要件

点列 $\{{\bf x}_n\}_{n=1}^N\ ({\bf x}_n\in{\mathbb R}^m)$ にフィットする3次ベジェ曲線を見つけたいとします。
点列の数 $N$ は4以上であるとします。(後で2以上でも導出可能なアイデアを下に書いております。)
3次ベジェ曲線の制御点を ${\bf P}_0,{\bf P}_1,{\bf P}_2,{\bf P}_3\in{\mathbb R}^m$ とします。
求める3次ベジェ曲線の端点 ${\bf P}_0,{\bf P}_3$ はそれぞれ点列の始点 ${\bf x}_1$ と終点 ${\bf x}_N$ に一致するとします。

$\begin{eqnarray} \left\{ \begin{array}{l} {\bf P}_0={\bf x}_1\\ {\bf P}_3={\bf x}_N \end{array} \right.\tag{1} \end{eqnarray}$

3次ベジェ曲線上の点を ${\bf P}(t)\ (t\in[0,1])$ とおくと以下のようになります。

$\begin{eqnarray} {\bf P}(t)=(1-t)^3{\bf P}_0+3(1-t)^2t{\bf P}_1+3(1-t)t^2{\bf P}_2+t^3{\bf P}_3\tag{2} \end{eqnarray}$

${\bf P}_0,{\bf P}_3$ はすでに求まっているので、 ${\bf P}_1,{\bf P}_2$ を求めれば、3次ベジェ曲線が定まることになります。

最小二乗法で解く

一般的に3次ベジェ曲線と点の距離を計算するのは大変なので、点に対応する3次ベジェ曲線のパラメータの値を計算します。
隣り合う点との距離の和の比を3次ベジェ曲線のパラメータの値とします。

$\begin{eqnarray} t_n= \left\{ \begin{array}{l} 0,\ (n = 1)\\ \dfrac{\displaystyle\sum_{i=1}^{n-1} || {\bf x}_{i+1}-{\bf x}_i ||}{\displaystyle\sum_{i=1}^{N-1} || {\bf x}_{i+1}-{\bf x}_i ||},\ (n \not=1,N)\\ 1.\ (n= N) \end{array} \right.\tag{3} \end{eqnarray}$

点列と3次ベジェ曲線との距離の2乗の和 ${\rm S}({\bf P}_1,{\bf P}_2)$ を計算します。

$\begin{eqnarray} {\rm S}({\bf P}_1,{\bf P}_2)&=&\sum_{n=1}^N || {\bf P}(t_n) - {\bf x}_n ||^2\\ &=&\sum_{n=1}^N{\rm L}_n\tag{4}\\ &=&\sum_{n=1}^N\sum_{i = 1}^m{\rm L}_{ni}\tag{5} \end{eqnarray}$

式 $(4)$ で ${\rm L}_n=|| {\bf P}(t_n) - {\bf x}_n ||^2$ とおきました。
ここで ${\bf P}_k,{\bf x}_n$ の第 $i$ 成分をそれぞれ ${\rm P}_{ki},x_{ni}$ と書くことにします。

式 $(5)$ で ${\rm L}_n$ はノルムの2乗なので、各成分の2乗の和となっています。
ここでわざわざ ${\rm L}_{ni}$ という関数の和にしたのは、この後に ${\rm P}_{1i},{\rm P}_{2i}$ で偏微分するのですが、
その時成分が異なる項は0になるためです。

$\begin{eqnarray} \dfrac{\partial {\rm L}_n}{\partial{\rm P}_{ki}}=\dfrac{\partial {\rm L}_{ni}}{\partial{\rm P}_{ki}}\tag{6}\\ \end{eqnarray}$

${\rm L}_{ni}$ を整理します。

$\begin{eqnarray} {\rm L}_{ni}&=&\left\{(1-t_n)^3{\rm P}_{0i} + 3(1-t_n)^2t_n{\rm P}_{1i} + 3(1-t_n)t_n^2{\rm P}_{2i} + t_n^3{\rm P}_{3i}-x_{ni}\right\}^2\\ &=&\left\{\underbrace{3(1-t_n)^2t_n}_{=a_n}{\rm P}_{1i} + \underbrace{3(1-t_n)t_n^2}_{=b_n}{\rm P}_{2i} + \underbrace{(1-t_n)^3{\rm P}_{0i} + t_n^3{\rm P}_{3i}-x_{ni}}_{=c_{ni}}\right\}^2\\ &=&\left(a_n{\rm P}_{1i}+b_n{\rm P}_{2i}+c_{ni}\right)^2\\ &=&a_n^2{\rm P}_{1i}^2+2a_nc_{ni}{\rm P}_{1i}+b_n^2{\rm P}_{2i}^2+2b_nc_{ni}{\rm P}_{2i}+2a_nb_n{\rm P}_{1i}{\rm P}_{2i} + c_{ni}^2\tag{7} \end{eqnarray}$

$\dfrac{\partial {\rm L}_{ni}}{\partial{\rm P}_{1i}}$ を計算します。

$\begin{eqnarray} \frac{\partial {\rm L}_{ni}}{\partial{\rm P}_{1i}}&=&\frac{\partial}{\partial{\rm P}_{1i}}\left(a_n^2{\rm P}_{1i}^2+2a_nc_{ni}{\rm P}_{1i}+b_n^2{\rm P}_{2i}^2+2b_nc_{ni}{\rm P}_{2i}+2a_nb_n{\rm P}_{1i}{\rm P}_{2i} + c_{ni}^2\right)\\ &=&2a_n^2{\rm P}_{1i}+2a_nb_n{\rm P}_{2i}+2a_nc_{ni}\tag{8} \end{eqnarray}$

$\dfrac{\partial {\rm L}_{ni}}{\partial{\rm P}_{2i}}$ を計算します。

$\begin{eqnarray} \frac{\partial {\rm L}_{ni}}{\partial{\rm P}_{2i}}&=&\frac{\partial}{\partial{\rm P}_{2i}}\left(a_n^2{\rm P}_{1i}^2+2a_nc_{ni}{\rm P}_{1i}+b_n^2{\rm P}_{2i}^2+2b_nc_{ni}{\rm P}_{2i}+2a_nb_n{\rm P}_{1i}{\rm P}_{2i} + c_{ni}^2\right)\\ &=&2a_nb_n{\rm P}_{1i}+2b_n^2{\rm P}_{2i}+2b_nc_{ni}\tag{9} \end{eqnarray}$

これでやっと ${\rm S}$ を偏微分する準備ができました。
解き方の基本方針としては、各成分ごとに方程式を解いて、値を求めていきます。

$\dfrac{\partial {\rm S}}{\partial{\rm P}_{1i}}=0$ を計算します。

$\begin{eqnarray} &&\frac{\partial {\rm S}}{\partial{\rm P}_{1i}}=0\\ &&\Rightarrow\sum_{n=1}^N\sum_{j = 1}^m\frac{\partial {\rm L}_{nj}}{\partial{\rm P}_{1i}}=0\\ &&\Rightarrow\sum_{n=1}^N\frac{\partial {\rm L}_{ni}}{\partial{\rm P}_{1i}}=0\tag{10} \end{eqnarray}$

式 $(8)$ を式 $(10)$ に代入します。

$\begin{eqnarray} &&\Rightarrow\sum_{n=1}^N\left(2a_n^2{\rm P}_{1i}+2a_nb_n{\rm P}_{2i}+2a_nc_{ni}\right)=0\\ &&\Rightarrow\left(\sum_{n=1}^Na_n^2\right){\rm P}_{1i}+\left(\sum_{n=1}^Na_nb_n\right){\rm P}_{2i}+\sum_{n=1}^Na_nc_{ni}=0\tag{11} \end{eqnarray}$

$\dfrac{\partial {\rm S}}{\partial{\rm P}_{2i}}=0$ を計算します。

$\begin{eqnarray} &&\frac{\partial {\rm S}}{\partial{\rm P}_{2i}}=0\\ &&\Rightarrow\sum_{n=1}^N\sum_{j = 1}^m\frac{\partial {\rm L}_{nj}}{\partial{\rm P}_{2i}}=0\\ &&\Rightarrow\sum_{n=1}^N\frac{\partial {\rm L}_{ni}}{\partial{\rm P}_{2i}}=0\tag{12} \end{eqnarray}$

式 $(9)$ を式 $(12)$ に代入します。

$\begin{eqnarray} &&\Rightarrow\sum_{n=1}^N\left(2a_nb_n{\rm P}_{1i}+2b_n^2{\rm P}_{2i}+2b_nc_{ni}\right)=0\\ &&\Rightarrow\left(\sum_{n=1}^Na_nb_n\right){\rm P}_{1i}+\left(\sum_{n=1}^Nb_n^2\right){\rm P}_{2i}+\sum_{n=1}^Nb_nc_{ni}=0\tag{13} \end{eqnarray}$

後は、式 $(11),(13)$ の連立方程式を解くだけです。

$\begin{eqnarray} \left\{ \begin{array}{l} \left(\displaystyle\sum_{n=1}^Na_n^2\right){\rm P}_{1i}+\left(\displaystyle\sum_{n=1}^Na_nb_n\right){\rm P}_{2i}+\displaystyle\sum_{n=1}^Na_nc_{ni}=0\\ \left(\displaystyle\sum_{n=1}^Na_nb_n\right){\rm P}_{1i}+\left(\displaystyle\sum_{n=1}^Nb_n^2\right){\rm P}_{2i}+\displaystyle\sum_{n=1}^Nb_nc_{ni}=0 \end{array} \right.\tag{14} \end{eqnarray}$

ここで、 $n=1,N$ のときは、それぞれ $t_n=0,1$ であり、このとき $a_n=0,b_n=0$ なので、式 $(14)$ の和は $2$ から $N-1$ まででよいことが分かります。

$\begin{eqnarray} \left\{ \begin{array}{l} \left(\displaystyle\sum_{n=2}^{N-1}a_n^2\right){\rm P}_{1i}+\left(\displaystyle\sum_{n=2}^{N-1}a_nb_n\right){\rm P}_{2i}+\displaystyle\sum_{n=2}^{N-1}a_nc_{ni}=0\\ \left(\displaystyle\sum_{n=2}^{N-1}a_nb_n\right){\rm P}_{1i}+\left(\displaystyle\sum_{n=2}^{N-1}b_n^2\right){\rm P}_{2i}+\displaystyle\sum_{n=2}^{N-1}b_nc_{ni}=0 \end{array} \right.\tag{15} \end{eqnarray}$

式 $(15)$ を行列を使って書きます。

$\begin{eqnarray} &&\begin{pmatrix} \displaystyle\sum_{n=2}^{N-1}a_n^2 & \displaystyle\sum_{n=2}^{N-1}a_nb_n\\ \displaystyle\sum_{n=2}^{N-1}a_nb_n & \displaystyle\sum_{n=2}^{N-1}b_n^2\\ \end{pmatrix} \begin{pmatrix} {\rm P}_{1i}\\ {\rm P}_{2i}\\ \end{pmatrix} = \begin{pmatrix} -\displaystyle\sum_{n=2}^{N-1}a_nc_{ni}\\ -\displaystyle\sum_{n=2}^{N-1}b_nc_{ni}\\ \end{pmatrix}\\ &&\Rightarrow \begin{pmatrix} {\rm P}_{1i}\\ {\rm P}_{2i}\\ \end{pmatrix} = \begin{pmatrix} \displaystyle\sum_{n=2}^{N-1}a_n^2 & \displaystyle\sum_{n=2}^{N-1}a_nb_n\\ \displaystyle\sum_{n=2}^{N-1}a_nb_n & \displaystyle\sum_{n=2}^{N-1}b_n^2\\ \end{pmatrix}^{-1} \begin{pmatrix} -\displaystyle\sum_{n=2}^{N-1}a_nc_{ni}\\ -\displaystyle\sum_{n=2}^{N-1}b_nc_{ni}\\ \end{pmatrix}\tag{16} \end{eqnarray}$

始点終点を除くデータが全て一致しなければ、式 $(16)$ の逆行列は存在します。(下でもう少し詳しく説明します。)

$\begin{eqnarray} \Rightarrow \begin{pmatrix} {\rm P}_{1i}\\ {\rm P}_{2i}\\ \end{pmatrix} &=& \frac{1}{\left(\displaystyle\sum_{n=2}^{N-1}a_n^2\right)\left(\displaystyle\sum_{n=2}^{N-1}b_n^2\right)-\left(\displaystyle\sum_{n=2}^{N-1}a_nb_n\right)^2} \begin{pmatrix} \displaystyle\sum_{n=2}^{N-1}b_n^2 & -\displaystyle\sum_{n=2}^{N-1}a_nb_n\\ -\displaystyle\sum_{n=2}^{N-1}a_nb_n & \displaystyle\sum_{n=2}^{N-1}a_n^2\\ \end{pmatrix} \begin{pmatrix} -\displaystyle\sum_{n=2}^{N-1}a_nc_{ni}\\ -\displaystyle\sum_{n=2}^{N-1}b_nc_{ni}\\ \end{pmatrix}\\ &=& \frac{1}{\left(\displaystyle\sum_{n=2}^{N-1}a_n^2\right)\left(\displaystyle\sum_{n=2}^{N-1}b_n^2\right)-\left(\displaystyle\sum_{n=2}^{N-1}a_nb_n\right)^2} \begin{pmatrix} -\left(\displaystyle\sum_{n=2}^{N-1}b_n^2\right)\left(\displaystyle\sum_{n=2}^{N-1}a_nc_{ni}\right) + \left(\displaystyle\sum_{n=2}^{N-1}a_nb_n\right)\left(\displaystyle\sum_{n=2}^{N-1}b_nc_{ni}\right)\\ \left(\displaystyle\sum_{n=2}^{N-1}a_nb_n\right)\left(\displaystyle\sum_{n=2}^{N-1}a_nc_{ni}\right) - \left(\displaystyle\sum_{n=2}^{N-1}a_n^2\right)\left(\displaystyle\sum_{n=2}^{N-1}b_nc_{ni}\right)\\ \end{pmatrix}\tag{17} \end{eqnarray}$

式 $(17)$ より、 ${\rm P}_{1i},{\rm P}_{2i}$ が求まりました。
これを全成分求めれば、 ${\bf P}_1,{\bf P}_2$ が求まります。

お疲れさまでした。

点列を増やす

このアルゴリズムで点列の数が少ない場合、点列にフィットしすぎて(過学習)あまりよくありません。
下図は、データ数が4つの場合。

また、データ点が3個の時もどうにか3次ベジェ曲線を求めたいです。

そこで、データを増やすことにします。
点列の間に新たに点を作成する(隣り合う点の中点)ことにより、点列を水増しします。
点列がN個なら2N-1個になります。

下図は、データ数4つを7つに増やした場合。上図と比べて、ややいい感じになっています。

※この点列を増やすロジックは私の好みもありますので、自身の好みに応じて変更してください。

式(16)の逆行列が存在する理由

式 $(16)$ の逆行列を求めるときに必要な行列式は、 $\left(\displaystyle\sum_{n=2}^{N-1}a_n^2\right)\left(\displaystyle\sum_{n=2}^{N-1}b_n^2\right)-\left(\displaystyle\sum_{n=2}^{N-1}a_nb_n\right)^2$ であり、コーシーシュワルツの不等式を適用すると、以下の式が成り立ちます。

$\begin{eqnarray} \left(\displaystyle\sum_{n=2}^{N-1}a_n^2\right)\left(\displaystyle\sum_{n=2}^{N-1}b_n^2\right)≥\left(\displaystyle\sum_{n=2}^{N-1}a_nb_n\right)^2\tag{18} \end{eqnarray}$

式 $(18)$ の等号が成り立たなければ、行列式が $\not=0$ となり、逆行列が存在するわけです。
式 $(18)$ の等号が成り立つのは、 $\dfrac{b_2}{a_2}=\dfrac{b_3}{a_3}=\cdots=\dfrac{b_{N-1}}{a_{N-1}}$ のときです。( $a_n\not=0\ (n\in{\mathbb N},1\leq n \leq N-1)$ )

ここで、 $\dfrac{b_n}{a_n}$ を計算してみます。

$\begin{eqnarray} \dfrac{b_n}{a_n}=(1-t_n)^2t_n/(1-t_n)t_n^2=(1-t_n)/t_n=(1-t_n)/t_n\tag{19} \end{eqnarray}$

横軸が $t_n$ で、縦軸が $(1-t_n)/t_n$ のグラフ。

$t_n$ は式 $(3)$ による定義より、同一点が連続しなければ $n>m\Rightarrow t_n>t_m$ です。
加えて、式 $(19)$ は単射なので、式 $(18)$ の等号が成り立つのは、始点と終点を除いたデータが全て一致するときのみであることが分かります。
そのような不正なデータは本アルゴリズムを適用する前に避ければよいので、特に問題にはなりません。

よって、式 $(16)$ の逆行列は(よっぽど変なデータでない限り)存在します。

成分ごとに解を求める理由

求めるベクトル ${\bf P}_1,{\bf P}_2$ で $\rm S$ をまとめて偏微分すれば、いいんじゃないのか？と思っている人がいるかもしれません。
はい、それでも勿論解が求まります。
実際に $m=2$ の場合、すなわち、 ${\bf x}_n,{\bf P}_k\in{\mathbb R}^2$ の場合で計算をしてみます。
ここで、 ${\bf Q}=({\bf P}_1^\top\ {\bf P}_2^\top)^\top=({\rm P}_{11}\ {\rm P}_{12}\ {\rm P}_{21}\ {\rm P}_{22})^\top\in{\mathbb R}^4$ とおきます。
$\dfrac{\partial{\rm S}}{\partial{\bf Q}}={\bf 0}$ と置きます。

$\begin{eqnarray} &&\dfrac{\partial{\rm S}}{\partial{\bf Q}}={\bf 0}\\ &&\vdots\\ &&\Rightarrow \begin{pmatrix} s_{11} & s_{12} & 0 & 0\\ s_{21} & s_{22} & 0 & 0\\ 0 & 0 & s_{33} & s_{34}\\ 0 & 0 & s_{43} & s_{44}\\ \end{pmatrix} \begin{pmatrix} {\rm P}_{11}\\ {\rm P}_{12}\\ {\rm P}_{21}\\ {\rm P}_{22}\\ \end{pmatrix} = \begin{pmatrix} {\rm v}_1\\ {\rm v}_2\\ {\rm v}_3\\ {\rm v}_4 \end{pmatrix}\tag{20}\\ &&\Rightarrow\left\{ \begin{array}{l} \begin{pmatrix} s_{11} & s_{12}\\ s_{21} & s_{22}\\ \end{pmatrix} \begin{pmatrix} {\rm P}_{11}\\ {\rm P}_{12} \end{pmatrix} =\begin{pmatrix} {\rm v}_1\\ {\rm v}_2 \end{pmatrix}\\ \begin{pmatrix} s_{33} & s_{34}\\ s_{43} & s_{44}\\ \end{pmatrix} \begin{pmatrix} {\rm P}_{21}\\ {\rm P}_{22} \end{pmatrix} =\begin{pmatrix} {\rm v}_3\\ {\rm v}_4 \end{pmatrix} \end{array} \right.\tag{21}\\ \end{eqnarray}$

式 $(20)$ の行列を見てもらえばわかる通り、無関係な成分の行列の要素が $0$ になっているので、式 $(19)$ のように変形できます。
ですので、各成分ごとに解を求めることにしたわけです。
式 $(20)$ の $s_{ij},{\rm v}_{i}$ は特に意味はありません。

少しでも早く動かす

$\displaystyle\sum_{n=2}^{N-1}a_n^2, \displaystyle\sum_{n=2}^{N-1}a_nb_n, \displaystyle\sum_{n=2}^{N-1}b_n^2$ は $t_n$ を求めた時点で計算可能なので、
事前に計算しておくと、少し計算時間が短くなります。

パラメータについて少し考察

$t_n$ を ${\mathbb R}^m$ のノルムで計算してますが、各成分ごとに $t_{ni}$ を求め、それを使って距離を計算するのもありかもしれません。

最後に

3次ベジェ曲線の制御点を ${\bf P}_0,{\bf P}_1,{\bf P}_2,{\bf P}_3$ と大文字の太字で書いていますが、(ベクトルは普通、小文字の太字)
私が参照しているベジェ曲線の文献では、ベクトルを大文字で書いているものが多く、それに倣ってその記法ですでに記事を書いているためです。

機械学習基礎理論独習

誤りがあればご指摘いただけると幸いです。数式が整うまで少し時間かかります。リンクフリーです。