2次元の点列にフィットする3次ベジェ曲線を求める

やりたいこと

点列 $\{(x_n,y_n)\in{\mathbb R}^2\}_{n=1}^N$ が与えられたとき、それにフィットする3次ベジェ曲線を1つ見つけたいとします。
要は、3次ベジェ曲線の4つの制御点の位置ベクトル $\{{\bf p}_i\in{\mathbb R}^2\}_{i = 0}^3$ を求めればよいわけです。

問題を書き換えてみる

3次ベジェ曲線の位置ベクトルを ${\bf r}(t)=(x(t)\ y(t))^\top$ とすると、 $x(t),y(t)$ は以下のようになります。
$t$ は3次ベジェ曲線のパラメータで $t\in[0,1]$ とします。

$\begin{eqnarray} &&x(t)=(1-t)^3p_{0x}+3(1-t)^2tp_{1x}+3(1-t)t^2p_{2x}+t^3p_{3x}\tag{1}\\ &&y(t)=(1-t)^3p_{0y}+3(1-t)^2tp_{1y}+3(1-t)t^2p_{2y}+t^3p_{3y}\tag{2}\\ \end{eqnarray}$

式 $(1),(2)$ は ${\bf p}_i$ の $x,y$ 成分をそれぞれ $p_{ix}, p_{iy}$ としました。
式 $(1),(2)$ を見てわかる通り $x$ と $y$ を入れ替えると式が同じなので、以降 $x(t)$ についてのみ考えることにします。
また、 $p_{ix}$ を $p_i$ と書くことにします。

式 $(1)$ を変形します。

$\begin{eqnarray} x(t)&=&(1-t)^3p_0+3(1-t)^2tp_1+3(1-t)t^2p_2+t^3p_3\\ &=& \begin{pmatrix}p_0 & p_1 & p_2 & p_3\end{pmatrix} \begin{pmatrix}(1-t)^3\\ 3(1-t)^2t\\ 3(1-t)t^2\\ t^3\end{pmatrix}\\ &=& {\bf p}^\top \begin{pmatrix}\phi_0(t)\\ \phi_1(t)\\ \phi_2(t)\\ \phi_3(t)\end{pmatrix}\\ &=&{\bf p}^\top{\boldsymbol\phi}(t)\tag{3} \end{eqnarray}$

式 $(3)$ の導出の過程で、 ${\bf p}:=\begin{pmatrix}p_0\\p_1\\p_2\\p_3\end{pmatrix},\ \phi_0(t):=(1-t)^3,\ \phi_1(t):=3(1-t)^2t,\ \phi_2(t):=3(1-t)t^2,\ \phi_3(t):=t^3, {\boldsymbol\phi}(t):=\begin{pmatrix}\phi_0(t)\\ \phi_1(t)\\ \phi_2(t)\\ \phi_3(t)\end{pmatrix}$ とおきました。
点列 $\{(t_n,x_n)\}_{n=1}^N$ を作成します。 $t_n$ は以下のように定義します。

$\begin{eqnarray} t_n= \left\{ \begin{array}{l} \dfrac{\displaystyle\sum_{i=2}^n\sqrt{(x_i-x_{i-1})^2+(y_i-y_{i-1})^2}}{\displaystyle\sum_{i=2}^N\sqrt{(x_i-x_{i-1})^2+(y_i-y_{i-1})^2}}\ \ \ (n > 1)\\ 0\hspace{ 300px } (n = 1) \end{array} \right.\tag{6} \end{eqnarray}$

式 $(6)$ で、 $x,y$ の両方の値を参照して $2$ 次元の距離を用いていることに注意してください。

以上より、問題は以下のように書き直すことができます。

「 $\{(t_n,x_n)\}_{n=1}^N$ が与えられたとき、それらに当てはまる曲線 $x(t)={\bf p}^\top{\boldsymbol\phi}(t)$ を求めよ。」

これは、いわゆる「線形回帰」の問題であることが分かります。

最小二乗法で解く

目的関数 ${\rm E}$ を定義します。

$\begin{eqnarray} {\rm E}({\bf p})&=&\sum_{n=1}^N\left(x(t_n)-x_n\right)^2\\ &=&\sum_{n=1}^N\left({\bf p}^\top{\boldsymbol\phi}(t_n)-x_n\right)^2\tag{7} \end{eqnarray}$

ここで、 ${\boldsymbol\Phi}=\begin{pmatrix}{\boldsymbol\phi}(t_1)^\top \\ {\boldsymbol\phi}(t_2)^\top \\ \vdots\\ {\boldsymbol\phi}(t_N)^\top \end{pmatrix},\ {\bf x}=\begin{pmatrix}x_1 \\ x_2 \\ \vdots \\ x_N\end{pmatrix}$ とおくと、式 $(7)$ は以下のように書き直せます。

$\begin{eqnarray} {\rm E}({\bf p})&=&||{\boldsymbol\Phi}{\bf p}-{\bf x}||^2\tag{8} \end{eqnarray}$

${\rm E}$ を ${\bf p}$ で微分して $={\bf 0}$ とおきます。

$\begin{eqnarray} &&\dfrac{\partial E}{\partial{\bf p}}={\bf 0}\\ &&\Rightarrow {\bf p}=({\boldsymbol\Phi}^\top{\boldsymbol\Phi})^{-1}{\boldsymbol\Phi}^\top{\bf x}\tag{9} \end{eqnarray}$

式 $(9)$ の導出については、曲線フィッティング - 最小二乗法の記事を参考にしてください。
正則化する場合は、リッジ回帰(L2正則化) の記事を参考にしてください。

式 $(9)$ で求めた $\bf p$ は、 ${\bf p}=\begin{pmatrix}p_{0x}\\p_{1x}\\p_{2x}\\p_{3x}\end{pmatrix}$ です。
$y$ についても同様に、点列 $\{(t_n,y_n)\}_{n=1}^N$ から式 $(9)$ を用いて ${\bf p}=\begin{pmatrix}p_{0y}\\p_{1y}\\p_{2y}\\p_{3y}\end{pmatrix}$ を求めます。

以上より、3次ベジェ曲線の4つの制御点の位置ベクトル $\{{\bf p}_i\in{\mathbb R}^2\}_{i = 0}^3$ が求まりました。

弧長パラメータの更新

式 $(6)$ を見てわかる通り、 $t_n$ は点列の弧長パラメータです。
これを曲線のパラメータとみなして、最小二乗法で3次ベジェ曲線を求めたわけです。
今度は、この求めた3次ベジェ曲線の幾何的な性質を用いて、 $t_n$ 更新します。

点 ${\bf x}_n$ から曲線へ下した垂線と曲線の交点が ${\bf r}(t_n)$ であるとき、以下の式が成り立ちます。

$\begin{eqnarray} ({\bf r}(t_n)-{\bf x}_n)^\top {\bf r}'(t_n)=0\tag{10} \end{eqnarray}$

式 $(10)$ を満たすような $t_n$ を求めるのですが、解析に解くのは困難なのでニュートン法で解きます。
式を簡単にするため、 $f(t)=({\bf r}(t)-{\bf x}_n)^\top {\bf r}'(t),\ {\bf r}_1(t)={\bf r}(t)-{\bf x}_n,\ {\bf r}_2={\bf r}'(t)$ とおくと、式 $(10)$ は以下のように書けます。

$\begin{eqnarray} f(t_n)={\bf r}_1(t_n)^\top {\bf r}_2(t_n)=0\tag{11} \end{eqnarray}$

式 $(11)$ をニュートン法で解くには、

$\begin{eqnarray} t_n\leftarrow t_n - \dfrac{f(t_n)}{f'(t_n)}\tag{12} \end{eqnarray}$

と $t_n$ を更新すればよいです。

$f'(t)=\dfrac{{\rm d}f}{{\rm d}t}$ を計算します。

$\begin{eqnarray} \dfrac{{\rm d}f}{{\rm d}t}&=&\dfrac{{\rm d}}{{\rm d}t} ({\bf r}(t)-{\bf x}_n)^\top {\bf r}'(t)\\ &=& \dfrac{{\rm d}}{{\rm d}t}\left({\bf r}(t)^\top {\bf r}'(t)-{\bf x}_n^\top {\bf r}'(t)\right)\\ &=& \dfrac{{\rm d}}{{\rm d}t}{\bf r}(t)^\top {\bf r}'(t)-\dfrac{{\rm d}}{{\rm d}t}{\bf x}_n^\top {\bf r}'(t)\\ &=& {\bf r}(t)'^\top {\bf r}'(t)+{\bf r}(t)^\top {\bf r}''(t)-{\bf x}_n^\top {\bf r}''(t)\\ &=& || {\bf r}'(t) ||^2+({\bf r}(t)-{\bf x}_n)^\top {\bf r}''(t)\tag{13} \end{eqnarray}$

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

$\begin{eqnarray} t_n\leftarrow t_n - \dfrac{({\bf r}(t_n)-{\bf x}_n)^\top {\bf r}'(t_n)}{ || {\bf r}'(t_n) ||^2+({\bf r}(t_n)-{\bf x}_n)^\top {\bf r}''(t_n)}\tag{14} \end{eqnarray}$

式 $(14)$ を用いて全ての $t_n\ (n=1,2,\ldots,N)$ を更新します。(ニュートン法は普通複数回行いますが、ここは1回でよいと思います。)
この時 $t_n\not\in[0,1]$ となる可能性があります。それは本記事では問題ありませんが、念のため記載しておきます。

その後、式 $(9)$ を用いて3次ベジェ曲線を更新するのですが、その時に誤差関数 $\rm E$ の値が必ず小さくなるわけではないので注意しましょう。
ですので、以下に示すような最大誤差や平均二乗平方根誤差で、3次ベジェ曲線を更新するかどうか決定したほうが良いと思います。

最大誤差と平均二乗平方根誤差

${\bf p}_x=\begin{pmatrix}p_{0x}\\p_{1x}\\p_{2x}\\p_{3x}\end{pmatrix}\in{\mathbb R}^4,\ {\bf p}_y=\begin{pmatrix}p_{0y}\\p_{1y}\\p_{2y}\\p_{3y}\end{pmatrix}\in{\mathbb R}^4$ とおくと、最大誤差 ${\rm E}_{\max}$ は以下のようになります。

$\begin{eqnarray} {\rm E}_{\max}&=&\max_{n\in\{1,\cdots,N\}}\sqrt{\left({\bf p}_x^\top{\boldsymbol\phi}(t_n)-x_n\right)^2+\left({\bf p}_y^\top{\boldsymbol\phi}(t_n)-y_n\right)^2}\tag{15} \end{eqnarray}$

${\bf P}=\begin{pmatrix}{\bf p}_x & {\bf p}_y\end{pmatrix}\in{\mathbb R}^{4\times 2},\ {\bf x}_n=\begin{pmatrix}x_n\\y_n\end{pmatrix}\in{\mathbb R}^2$ とおくと、式 $(15)$ は以下のように書けます。

$\begin{eqnarray} {\rm E}_{\max}&=&\max_{n\in\{1,\cdots,N\}}||{\bf P}^\top{\boldsymbol\phi}(t_n)-{\bf x}_n||\tag{16} \end{eqnarray}$

最大誤差 ${\rm E}_{\max}$ の大小のみを比較するのなら、式 $(15)$ の平方根は不要ですし、式 $(16)$ のノルムは2乗してもよいと思います。

平均二乗平方根誤差 ${\rm E}_{\rm RMS}$ とは、「誤差の二乗の平均の平方根」の事なので以下のようになります。

$\begin{eqnarray} {\rm E}_{\rm RMS}&=&\sqrt{\dfrac{1}{N}\left({\rm E}({\bf p}_x)+{\rm E}({\bf p}_y)\right)}\\ &=&\dfrac{1}{\sqrt{N}}\sqrt{\sum_{n=1}^N\left({\bf p}_x^\top{\boldsymbol\phi}(t_n)-x_n\right)^2+\sum_{n=1}^N\left({\bf p}_y^\top{\boldsymbol\phi}(t_n)-y_n\right)^2}\\ &=&\dfrac{1}{\sqrt{N}}\sqrt{\sum_{n=1}^N\left(\left({\bf p}_x^\top{\boldsymbol\phi}(t_n)-x_n\right)^2+\left({\bf p}_y^\top{\boldsymbol\phi}(t_n)-y_n\right)^2\right)}\\ &=&\dfrac{1}{\sqrt{N}}\sqrt{\sum_{n=1}^N||{\bf P}^\top{\boldsymbol\phi}(t_n)-{\bf x}_n||^2}\tag{17} \end{eqnarray}$

${\bf X}=\begin{pmatrix}{\bf x}_1^\top\\ {\bf x}_2^\top\\ \vdots\\ {\bf x}_N\end{pmatrix}$ とおくと、式 $(17)$ は以下のように書けます。

$\begin{eqnarray} {\rm E}_{\rm RMS}&=&\dfrac{1}{\sqrt{N}}||{\boldsymbol\Phi}{\bf P}-{\bf X}||\tag{18} \end{eqnarray}$

平均二乗平方根誤差 ${\rm E}_{\rm RMS}$ の大小のみを比較するのなら、式 $(17)$ の $\dfrac{1}{\sqrt{N}}$ や平方根は不要ですし、式 $(18)$ の $\dfrac{1}{\sqrt{N}}$ も不要でノルムは2乗してもよいと思います。

まとめ

①：式 $(6)$ で $t_n$ を求める
②：式 $(9)$ で 3次ベジェ曲線を求める
以下③、④を2～4回繰り返す。
③：式 $(14)$ で $t_n$ を更新する
④：式 $(9)$ で3次ベジェ曲線を更新する
④完了時に、最大誤差 ${\rm E}_{\max}$ や平均二乗平方根誤差 ${\rm E}_{\rm RMS}$ の値が増えるようなら④の前の3次ベジェ曲線を採用し、処理を終了する。

参考リンク

曲線フィッティング - 最小二乗法
 リッジ回帰(L2正則化)

機械学習基礎理論独習

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