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

要件

2つの点列 $\{{\bf x}_m\in{\mathbb R}^2\}_{m=1}^{\rm M},\ \{{\bf y}_n\in{\mathbb R}^2\}_{n=1}^{\rm N}\ ({\rm M},{\rm N}\in{\mathbb N})$ が与えられていて、点列の終点と始点が一致( ${\bf x}_{\rm M}={\bf y}_1$ )します。
この2つの点列にフィットする2つのG1連続の3次ベジェ曲線 ${\bf p}:{\mathbb R}\rightarrow{\mathbb R}^2,\ {\bf q}:{\mathbb R}\rightarrow{\mathbb R}^2$ を求めたいとします。
${\bf p},{\bf q}$ の制御点をそれぞれ ${\bf p}_i,{\bf q}_i\ (i=0,1,2,3)$ とします。
G1連続なので ${\bf p}'(1)$ と ${\bf q}'(0)$ の方向が逆向き、すなわち ${\bf p}_2,{\bf p}_3,{\bf q}_1$ が一直線上に順番に並んでいる必要があります。

$\begin{aligned} \left\{ \begin{array}{l} {\bf q}_0={\bf p}_3\\ {\bf q}_1={\bf p}_3+ l\dfrac{{\bf p}_3-{\bf p}_2}{||{\bf p}_3-{\bf p}_2||}\\ l\geq 0\\ \end{array} \right.\ \cdots\ (1) \end{aligned}$

また、 ${\bf p}$ の始点は ${\bf x}_1$ と、 ${\bf q}$ の終点は ${\bf y}_{\rm N}$ と一致するものとします。

$\begin{aligned} \left\{ \begin{array}{l} {\bf p}(0)={\bf p}_0={\bf x}_1\\ {\bf q}(1)={\bf q}_3={\bf y}_{\rm N}\\ \end{array} \right.\ \cdots\ (2) \end{aligned}$

パラメータをまとめて ${\boldsymbol\theta}_p:=\{{\bf p}_1,{\bf p}_2,{\bf p}_3\},\ {\boldsymbol\theta}_q:=\{{\bf p}_2,{\bf p}_3,l,{\bf q}_2\},\ {\boldsymbol\theta}:={\boldsymbol\theta}_p\cup{\boldsymbol\theta}_q=\{{\bf p}_1,{\bf p}_2,{\bf p}_3,l,{\bf q}_2\}$ とおく。

${\bf p}({\boldsymbol\theta}_p;{\bf p}_0,t)$ の定義

$\begin{aligned} {\bf p}({\boldsymbol\theta}_p;{\bf p}_0,t) &={\bf p}({\bf p}_1,{\bf p}_2,{\bf p}_3;{\bf p}_0,t)\\ &=\displaystyle\sum_{i=0}^3{\rm B}_i^3(t){\bf p}_i\\ &={\rm B}_0^3(t){\bf p}_0+{\rm B}_1^3(t){\bf p}_1+{\rm B}_2^3(t){\bf p}_2+{\rm B}_3^3(t){\bf p}_3\ \cdots\ (3)\\ \end{aligned}$

${\bf q}({\boldsymbol\theta}_q;{\bf q}_3,u)$ の定義

$\begin{aligned} {\bf q}({\boldsymbol\theta}_q;{\bf q}_3,u) &={\bf q}({\bf p}_2,{\bf p}_3, l,{\bf q}_2;{\bf q}_3,u)\\ &=\displaystyle\sum_{i=0}^3{\rm B}_i^3(u){\bf q}_i\\ &={\rm B}_0^3(u){\bf q}_0+{\rm B}_1^3(u){\bf q}_1+{\rm B}_2^3(u){\bf q}_2+{\rm B}_3^3(u){\bf q}_3\\ &={\rm B}_0^3(u){\bf p}_3+{\rm B}_1^3(u)\left({\bf p}_3+ l\dfrac{{\bf p}_3-{\bf p}_2}{||{\bf p}_3-{\bf p}_2||}\right)+{\rm B}_2^3(u){\bf q}_2+{\rm B}_3^3(u){\bf q}_3\\ &=\dfrac{{\rm B}_1^3(u)l}{||{\bf p}_3-{\bf p}_2||}\left({\bf p}_3-{\bf p}_2\right)+\left({\rm B}_0^3(u)+{\rm B}_1^3(u)\right){\bf p}_3+{\rm B}_2^3(u){\bf q}_2+{\rm B}_3^3(u){\bf q}_3\ \cdots\ (4)\\ \end{aligned}$

要件を整理します。
「 $\{{\bf x}_m\}_{m=1}^{\rm M},\ \{{\bf y}_n\}_{n=1}^{\rm N}$ に当てはまるようなG1連続の ${\bf p}(t),\ {\bf q}(s)$ を求めるために ${\bf p}_1,{\bf p}_2,{\bf p}_3, l,{\bf q}_2$ を最適化せよ。」となります。

弧長パラメータの初期化

点列 $\{{\bf x}_m\}_{m=1}^{\rm M},\ \{{\bf y}_n\}_{n=1}^{\rm N}$ から弧長パラメータ $t_m\in[0,1],u_n\in[0,1]$ を生成します。
弧長パラメータを初期化についてはこちらの記事をご覧ください。

目的関数

目的関数 ${\rm E}({\boldsymbol\theta};{\bf p}_0,{\bf q}_3)$ の定義

$\begin{aligned} {\rm E}({\boldsymbol\theta};{\bf p}_0,{\bf q}_3) &={\rm E}_p({\boldsymbol\theta}_p;{\bf p}_0)+{\rm E}_q({\boldsymbol\theta}_q;{\bf q}_3)\\ &=\dfrac{1}{2}\displaystyle\sum_{m=2}^{{\rm M}}||{\bf p}({\boldsymbol\theta}_p;{\bf p}_0,t_m)-{\bf x}_m||^2+\dfrac{1}{2}\displaystyle\sum_{n=2}^{{\rm N}-1}||{\bf q}({\boldsymbol\theta}_q;{\bf q}_3,u_n)-{\bf y}_n||^2\ \cdots\ (5) \end{aligned}$

Eの微分

${\rm E}$ に関するもろもろの計算を先に済ませておきます。

$\dfrac{\partial{\rm E}_p}{\partial{\bf p}_i}\ (i\in\{1,2,3\})$ の計算

$\begin{aligned} \dfrac{\partial{\rm E}_p}{\partial{\bf p}_i} &=\dfrac{1}{2}\displaystyle\sum_{m=2}^{{\rm M}}\dfrac{\partial}{\partial{\bf p}_i}||{\bf p}({\boldsymbol\theta}_p;{\bf p}_0,t_m)-{\bf x}_m||^2\\ &=\dfrac{1}{2}\displaystyle\sum_{m=2}^{{\rm M}}\dfrac{\partial}{\partial{\bf p}_i}||{\rm B}_1^3(t_m){\bf p}_1+\cdots-{\bf x}_m||^2\\ &=\displaystyle\sum_{m=2}^{{\rm M}}{\rm B}_i^3(t_m)\left({\bf p}({\boldsymbol\theta}_p;{\bf p}_0,t_m)-{\bf x}_m\right)\ \cdots\ (6) \end{aligned}$

$\dfrac{\partial{\rm E}_q}{\partial{\bf q}_2}$ の計算

$\begin{aligned} \dfrac{\partial{\rm E}_q}{\partial{\bf q}_2} &=\displaystyle\sum_{n=2}^{{\rm N}-1}{\rm B}_2^3(u_n)\left({\bf q}({\bf p}_2,{\bf p}_3, l,{\bf q}_2;{\bf q}_3,u_n)-{\bf y}_n\right)\ \cdots\ (7) \end{aligned}$

$\dfrac{\partial{\rm E}_q}{\partial l}$ の計算

$\begin{aligned} \dfrac{\partial{\rm E}_q}{\partial l} &=\dfrac{1}{2}\displaystyle\sum_{n=2}^{{\rm N}-1}\dfrac{\partial}{\partial l}||{\bf q}({\boldsymbol\theta}_q;{\bf q}_3,u_n)-{\bf y}_n||^2\\ &=\dfrac{1}{2}\displaystyle\sum_{n=2}^{{\rm N}-1}\dfrac{\partial}{\partial l}\left|\left|{\rm B}_1^3(u_n) l\dfrac{{\bf p}_3-{\bf p}_2}{||{\bf p}_3-{\bf p}_2||}+\cdots-{\bf y}_n\right|\right|^2\\ &=\displaystyle\sum_{n=2}^{{\rm N-1}}{\rm B}_1^3(u_n)\dfrac{\left({\bf p}_3-{\bf p}_2\right)}{||{\bf p}_3-{\bf p}_2||}^\top\left({\bf q}({\boldsymbol\theta}_q;{\bf q}_3,u_n)-{\bf y}_n\right)\\ &=\dfrac{\left({\bf p}_3-{\bf p}_2\right)}{||{\bf p}_3-{\bf p}_2||}^\top\displaystyle\sum_{n=2}^{{\rm N-1}}{\rm B}_1^3(u_n)\left({\bf q}({\boldsymbol\theta}_q;{\bf q}_3,u_n)-{\bf y}_n\right)\ \cdots\ (8) \end{aligned}$

$\dfrac{\partial{\rm E}_q}{\partial{\bf p}_2}$ の計算

$\begin{aligned} \dfrac{\partial{\rm E}_q}{\partial{\bf p}_2} &=\dfrac{1}{2}\displaystyle\sum_{n=2}^{{\rm N}-1}\dfrac{\partial}{\partial{\bf p}_2}||{\bf q}({\boldsymbol\theta}_q;{\bf q}_3,u_n)-{\bf y}_n||^2\\ &=\dfrac{1}{2}\displaystyle\sum_{n=2}^{{\rm N}-1}\dfrac{\partial}{\partial{\bf p}_2}\left|\left|\dfrac{{\rm B}_1^3(u_n)l}{||{\bf p}_3-{\bf p}_2||}\left({\bf p}_3-{\bf p}_2\right)+\cdots-{\bf y}_n\right|\right|^2\\ &=\dfrac{1}{2}\displaystyle\sum_{n=2}^{{\rm N}-1}{\bf f}\big(\ \underbrace{{\rm B}_1^3(u_n)l,\ {\bf p}_3,\ {\bf p}_2,\ {\rm B}_0^3(u_n){\bf p}_3+{\rm B}_1^3(u_n){\bf p}_3+{\rm B}_2^3(u_n){\bf q}_2+{\rm B}_3^3(u_n){\bf q}_3-{\bf y}_n}_{=:{\boldsymbol\phi}_n}\ \big)\\ &=\dfrac{1}{2}\displaystyle\sum_{n=2}^{{\rm N}-1}{\bf f}\left({\boldsymbol\phi}_n\right)\ \cdots\ (9) \end{aligned}$

式 $(9)$ で ${\boldsymbol\phi}_n:=\{{\rm B}_1^3(u_n)l,\ {\bf p}_3,\ {\bf p}_2,\ {\rm B}_0^3(u_n){\bf p}_3+{\rm B}_1^3(u_n){\bf p}_3+{\rm B}_2^3(u_n){\bf q}_2+{\rm B}_3^3(u_n){\bf q}_3-{\bf y}_n\}$ とおきました。

式 $(9)$ で呼び出している ${\bf f}$ の定義

Matrix Calculus で norm2(a/norm2(x-y)*(x-y)+z)^2 を y で微分した結果

$\begin{aligned} &{\bf f}:\ {\mathbb R}^{1+2+2+2}\rightarrow{\mathbb R}^2,\\ &(a,{\bf x},{\bf y},{\bf z})\mapsto\dfrac{\partial}{\partial{\bf y}}\left|\left|\dfrac{a}{||{\bf x}-{\bf y}||}\left({\bf x}-{\bf y}\right)+{\bf z}\right|\right|^2=\dfrac{t_2}{t_1^3}{\bf t}_0^\top{\bf t}_3{\bf t}_0-\dfrac{t_2}{t_1}{\bf t}_3\\ &{\rm where}\ {\bf t}_0={\bf x}-{\bf y},\ t_1=||{\bf t}_0||,\ t_2=2a,\ {\bf t}_3={\bf z}+\dfrac{a}{t_1}{\bf t}_0 \end{aligned}$

$\dfrac{\partial{\rm E}_q}{\partial{\bf p}_3}$ の計算

$\begin{aligned} \dfrac{\partial{\rm E}_q}{\partial{\bf p}_3} &=\dfrac{1}{2}\displaystyle\sum_{n=2}^{{\rm N}-1}\dfrac{\partial}{\partial{\bf p}_3}||{\bf q}({\boldsymbol\theta}_q;{\bf q}_3,u_n)-{\bf y}_n||^2\\ &=\dfrac{1}{2}\displaystyle\sum_{n=2}^{{\rm N}-1}\dfrac{\partial}{\partial{\bf p}_3}\left|\left|\dfrac{{\rm B}_1^3(u_n)l}{||{\bf p}_3-{\bf p}_2||}\left({\bf p}_3-{\bf p}_2\right)+\left({\rm B}_0^3(u_n)+{\rm B}_1^3(u_n)\right){\bf p}_3+\cdots-{\bf y}_n\right|\right|^2\\ &=\dfrac{1}{2}\displaystyle\sum_{n=2}^{{\rm N}-1}{\bf g}\big(\ \underbrace{{\rm B}_1^3(u_n)l,\ {\rm B}_0^3(u_n)+{\rm B}_1^3(u_n),\ {\bf p}_3,\ {\bf p}_2,\ {\rm B}_2^3(u_n){\bf q}_2+{\rm B}_3^3(u_n){\bf q}_3-{\bf y}_n}_{=:{\boldsymbol\psi}_n}\ \big)\\ &=\dfrac{1}{2}\displaystyle\sum_{n=2}^{{\rm N}-1}{\bf g}\left({\boldsymbol\psi}_n\right)\ \cdots\ (10) \end{aligned}$

式 $(10)$ で ${\boldsymbol\psi}_n:=\{{\rm B}_1^3(u_n)l,\ {\rm B}_0^3(u_n)+{\rm B}_1^3(u_n),\ {\bf p}_3,\ {\bf p}_2,\ {\rm B}_2^3(u_n){\bf q}_2+{\rm B}_3^3(u_n){\bf q}_3-{\bf y}_n\}$ とおきました

式 $(10)$ で呼び出している ${\bf g}$ の定義

Matrix Calculus で norm2(a/norm2(x-y)*(x-y)+b*x+z)^2 を x で微分した結果

$\begin{aligned} &{\bf g}:\ {\mathbb R}^{1+1+2+2+2}\rightarrow{\mathbb R}^2,\\ &(a,b,{\bf x},{\bf y},{\bf z})\mapsto\dfrac{\partial}{\partial{\bf x}}\left|\left|\dfrac{a}{||{\bf x}-{\bf y}||}\left({\bf x}-{\bf y}\right)+b{\bf x}+{\bf z}\right|\right|^2=\dfrac{t_2}{t_1}{\bf t}_3-\dfrac{t_2}{t_1^3}{\bf t}_0^\top{\bf t}_3{\bf t}_0+2b{\bf t}_3\\ &{\rm where}\ {\bf t}_0={\bf x}-{\bf y},\ t_1=||{\bf t}_0||,\ t_2=2a,\ {\bf t}_3={\bf z}+\dfrac{a}{t_1}{\bf t}_0+b{\bf x} \end{aligned}$

勾配ベクトル

$\nabla{\rm E}$ の計算

$\begin{aligned} \nabla{\rm E} &=\dfrac{\partial{\rm E}}{\partial{\boldsymbol\theta}} =\begin{pmatrix} \dfrac{\partial{\rm E}}{\partial{\bf p}_1}\\ \dfrac{\partial{\rm E}}{\partial{\bf p}_2}\\ \dfrac{\partial{\rm E}}{\partial{\bf p}_3}\\ \dfrac{\partial{\rm E}}{\partial{\bf q}_2} \end{pmatrix} = \begin{pmatrix} \dfrac{\partial{\rm E}_p}{\partial{\bf p}_1}\\ \dfrac{\partial{\rm E}_p}{\partial{\bf p}_2}+\dfrac{\partial{\rm E}_q}{\partial{\bf p}_2}\\ \dfrac{\partial{\rm E}_p}{\partial{\bf p}_3}+\dfrac{\partial{\rm E}_q}{\partial{\bf p}_3}\\ \dfrac{\partial{\rm E}_q}{\partial{\bf q}_2} \end{pmatrix}\\ &= \begin{pmatrix} \displaystyle\sum_{m=2}^{{\rm M}}{\rm B}_1^3(t_m)\left({\bf p}({\boldsymbol\theta}_p;{\bf p}_0,t_m)-{\bf x}_m\right)\\ \displaystyle\sum_{m=2}^{{\rm M}}{\rm B}_2^3(t_m)\left({\bf p}({\boldsymbol\theta}_p;{\bf p}_0,t_m)-{\bf x}_m\right)+\dfrac{1}{2}\displaystyle\sum_{n=2}^{{\rm N}-1}{\bf f}\left({\boldsymbol\phi}_n\right)\\ \displaystyle\sum_{m=2}^{{\rm M}}{\rm B}_3^3(t_m)\left({\bf p}({\boldsymbol\theta}_p;{\bf p}_0,t_m)-{\bf x}_m\right)+\dfrac{1}{2}\displaystyle\sum_{n=2}^{{\rm N}-1}{\bf g}\left({\boldsymbol\psi}_n\right)\\ \dfrac{\left({\bf p}_3-{\bf p}_2\right)}{||{\bf p}_3-{\bf p}_2||}^\top\displaystyle\sum_{n=2}^{{\rm N-1}}{\rm B}_1^3(u_n)\left({\bf q}({\boldsymbol\theta}_q;{\bf q}_3,u_n)-{\bf y}_n\right)\\ \displaystyle\sum_{n=2}^{{\rm N}-1}{\rm B}_2^3(u_n)\left({\bf q}({\boldsymbol\theta}_q;{\bf q}_3,u_n)-{\bf y}_n\right) \end{pmatrix}\ \cdots\ (11) \end{aligned}$

弧長パラメータの更新

弧長パラメータの初期化⇒最適化⇒(弧長パラメータの更新⇒最適化)x2～4でフィットする2つのベジェ曲線が得られます。
弧長パラメータの更新についてはこちらの記事をご覧ください。

最後に

まだプログラムを組んでいないので、組んだらまた追記します。

機械学習基礎理論独習

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