はじめに

Graphics Gem1 の「AN ALGORITHM FOR AUTOMATICALLY FITTING DIGITIZED CURVES」(以下、元記事)を意訳していきます。
自分の意見は青色の太字で、書いていきます。
本記事でいうソースコードとは本アルゴリズムを実装した FitCurve.c ファイルの事です。
記号の書き方は元記事に合わせます。ベクトルが小文字の太字になってなかったりしますが、それは元記事に合わせているためです。
式番号は元論文と合わせていません。
なお、ベクトルの内積は ${\bf a}^\top{\bf b}$ と書くことにします。

登場する変数、定数

点列: $\{d_i\in{\mathbb R}^2\}_{i=1}^n$
指定誤差: $\epsilon(>0)$
粗めの誤差: $\Psi(>\epsilon)$
パラメータが紐づいた点列: $\{(u_i\in{\mathbb R},d_i\in{\mathbb R}^2)\}_{i=1}^n$
3次ベジェ曲線の位置ベクトル: $Q(t)\in{\mathbb R}^2\ (t\in[0,1])$
3次ベジェ曲線の制御点: $V_0,V_1,V_2,V_3\in{\mathbb R}^2$
$n$ 次のバーンスタイン基底関数: $B_i^n(t)=\begin{pmatrix}n\\ i\end{pmatrix}t^i(1-t)^{n-1}\hspace{30px}(i=0,\ldots,n)$
$V_0$ における単位接ベクトル: $\hat{t}_1\in{\mathbb R}^2$
$V_3$ における単位接ベクトル: $\hat{t}_2\in{\mathbb R}^2$
$V_0$ から $V_1$ までの距離: $\alpha_1(>0)$
$V_3$ から $V_2$ までの距離: $\alpha_2(>0)$

要件

「点列 $\{d_i\in{\mathbb R}^2\}_{i=1}^n$ が与えられたときに、指定した誤差 $\epsilon(>0)$ 以内の3次ベジェ曲線群を求める」というのが要件です。

①：3次ベジェ曲線の始点と終点の接ベクトルを計算する

3次ベジェ曲線の位置ベクトルを $Q(t)\in{\mathbb R}^2\ (t\in[0,1])$、制御点を $V_0,V_1,V_2,V_3\in{\mathbb R}^2$ とします。
$n$ 次のバーンスタイン基底関数を $B_i^n(t)=\begin{pmatrix}n\\i\end{pmatrix}t^i(1-t)^{n-1}\hspace{30px}(i=0,\ldots,n)$ と定めると、$Q(t)$ は以下のように書けます。

\begin{eqnarray}
Q(t)=\displaystyle\sum_{i=0}^3V_iB_i^3(t),\hspace{30px}t\in[0, 1]\tag{1}
\end{eqnarray}

$V_0,V_3$ の単位接ベクトルをそれぞれ $\hat{t_1},\hat{t_2}\in{\mathbb R}^2$ とします。
点列 $\{d_i\}_{i=1}^n$ から $\hat{t_1},\hat{t_2}$ を次のように定めます。

\begin{eqnarray}
\left\{
\begin{array}{l}
\hat{t}_1=\dfrac{d_1-d_0}{||d_1-d_0||}\\
\hat{t}_2=\dfrac{d_{n-1}-d_n}{||d_{n-1}-d_n||}
\end{array}
\right.\tag{2}
\end{eqnarray}

$\hat{t_1}$ は $\overrightarrow{V_0V_1}$ の方向を向いており、$\hat{t_2}$ は $\overrightarrow{V_3V_2}$ の方向を向いていることに注意しましょう。
この接ベクトルの決め方はやや荒いと思います。
元論文にもありますが、隣接点を1つ参照するのではなく、複数の隣接点までの平均ベクトルまたは
複数の隣接点から最小二乗法でフィットする直線を求め、その直線の方向ベクトルを
接ベクトルとして採用するとよいと思います。

②：点列に対応する曲線のパラメータの位置を計算する

後で誤差を計算するために、点列: $\{d_i\in{\mathbb R}^2\}_{i=1}^n$ に対応するパラメータ $u_i\in[0,1]$ を計算します。
このパラメータの求め方は粗いものであることを理解しておいてください。

\begin{eqnarray}
u_i=\left\{
\begin{array}{l}
\dfrac{\displaystyle\sum_{j=2}^i||d_{j}-d_{j-1}||}{\displaystyle\sum_{j=2}^n||d_{j}-d_{j-1}||}\hspace{30px}(i>1)\\
0\hspace{160px}(i=1)
\end{array}
\right.\tag{3}
\end{eqnarray}

ここの処理ですが、点列が $\{d_i\}_{i=1}^n$ から $\{(u_i,d_i)\}_{i=1}^n$ になるイメージです。

③：3次ベジェ曲線の始点と終点を定める

3次ベジェ曲線の始点(開始側の制御点) $V_0$ は点列の始点 $d_1$ と一致させます。
3次ベジェ曲線の終点(終了側の制御点) $V_3$ は点列の終点 $d_n$と一致させます。

\begin{eqnarray}
\left\{
\begin{array}{l}
V_0=d_1\\
V_3=d_n
\end{array}
\right.\tag{4}
\end{eqnarray}

④：始点と終点以外の2つ制御点の端点からの距離を求めて、3次ベジェ曲線を決定する

$V_0$ から $V_1$ までの距離を $\alpha_1(>0)$ 、$V_3$ から $V_2$ までの距離を $\alpha_2(>0)$ とおきます。
このとき、 $V_1,V_2$ は以下のように書けます。

\begin{eqnarray}
\left\{
\begin{array}{l}
V_1=\alpha_1\hat{t}_1+V_0\\
V_2=\alpha_2\hat{t}_2+V_3
\end{array}
\right.\tag{5}
\end{eqnarray}

3次ベジェ曲線 $Q$と点列 $\{d_i\}_{i=1}^n$ との距離の二乗和を $S$ とおきます。

\begin{eqnarray}
S&=&\sum_{i = 1}^n||d_i-Q(u_i)||^2\\
&=&\sum_{i = 1}^n(d_i-Q(u_i))^\top(d_i-Q(u_i))\tag{6}
\end{eqnarray}

$S$ を最小にするような $\alpha_1,\alpha_2$ を求めるために、$S$ を $\alpha_1,\alpha_2$ で偏微分して $=0$ とおきます。

\begin{eqnarray}
\left\{
\begin{array}{l}
\dfrac{\partial S}{\partial\alpha_1}=0\hspace{100px}(7)\\
\dfrac{\partial S}{\partial\alpha_2}=0\hspace{100px}(8)
\end{array}
\right.
\end{eqnarray}

式 $(7)$ の左辺を変形します。

\begin{eqnarray}
\dfrac{\partial S}{\partial\alpha_1}&=&\dfrac{\partial}{\partial\alpha_1}\sum_{i = 1}^n||d_i-Q(u_i)||^2\\
&=&\sum_{i = 1}^n2(d_i-Q(u_i))^\top\dfrac{\partial}{\partial\alpha_1}(d_i-Q(u_i))\\
&=&-2\sum_{i = 1}^n(d_i-Q(u_i))^\top\dfrac{\partial}{\partial\alpha_1}Q(u_i)\tag{9}
\end{eqnarray}

式 $(9)$ の $\dfrac{\partial}{\partial\alpha_1}Q(u_i)$ を計算します。

\begin{eqnarray}
\dfrac{\partial}{\partial\alpha_1}Q(u_i)&=&\dfrac{\partial}{\partial\alpha_1}\sum_{j=0}^3V_jB_j^3(u_i)\\
&=&\dfrac{\partial}{\partial\alpha_1}\left(V_0B_0^3(u_i)+V_1B_1^3(u_i)+V_2B_2^3(u_i)+V_3B_3^3(u_i)\right)\\
&=&\dfrac{\partial}{\partial\alpha_1}\left(V_0B_0^3(u_i)+\left(\alpha_1\hat{t}_1+V_0\right)B_1^3(u_i)+\left(\alpha_2\hat{t}_2+V_3\right)B_2^3(u_i)+V_3B_3^3(u_i)\right)\\
&=&\dfrac{\partial}{\partial\alpha_1}\alpha_1\hat{t}_1B_1^3(u_i)\\
&=&\hat{t}_1B_1^3(u_i)\tag{10}
\end{eqnarray}

式 $(9),(10)$ を式 $(7)$ に代入します。

\begin{eqnarray}
&&\dfrac{\partial S}{\partial\alpha_1}=0\\
&&\Rightarrow-2\sum_{i = 1}^n(d_i-Q(u_i))^\top\hat{t}_1B_1^3(u_i)=0\\
&&\Rightarrow\sum_{i = 1}^n(d_i-Q(u_i))^\top\hat{t}_1B_1^3(u_i)=0\\
&&\Rightarrow\sum_{i = 1}^n\underbrace{B_1^3(u_i)\hat{t}_1}_{:=A_{i,1}}^\top Q(u_i)=\sum_{i = 1}^n\underbrace{B_1^3(u_i)\hat{t}_1}_{:=A_{i,1}}^\top d_i\\
&&\Rightarrow\sum_{i = 1}^nA_{i,1}^\top Q(u_i)=\sum_{i = 1}^nA_{i,1}^\top d_i\tag{11}
\end{eqnarray}

式 $(11)$ で $A_{i,j}=B_j^3(u_i)\hat{t}_j$ とおきました。

式 $(11)$ の左辺を計算します。

\begin{eqnarray}
\sum_{i = 1}^nA_{i,1}^\top Q(u_i)&=&\sum_{i = 1}^nA_{i,1}^\top\left(V_0B_0^3(u_i)+\left(\alpha_1\hat{t}_1+V_0\right)B_1^3(u_i)+\left(\alpha_2\hat{t}_2+V_3\right)B_2^3(u_i)+V_3B_3^3(u_i)\right)\\
&=&\sum_{i = 1}^nA_{i,1}^\top\big(V_0B_0^3(u_i)+\alpha_1\underbrace{\hat{t}_1B_1^3(u_i)}_{=A_{i,1}}+V_0B_1^3(u_i)+\alpha_2\underbrace{\hat{t}_2B_2^3(u_i)}_{=A_{i,2}}+V_3B_2^3(u_i)+V_3B_3^3(u_i)\big)\\
&=&\underbrace{\left(\sum_{i = 1}^n||A_{i,1}||^2\right)}_{=:c_{1,1}}\alpha_1+\underbrace{\left(\sum_{i = 1}^nA_{i,1}^\top A_{i,2}\right)}_{=:c_{1,2}}\alpha_2+\sum_{i = 1}^nA_{i,1}^\top\big(V_0B_0^3(u_i)+V_0B_1^3(u_i)+V_3B_2^3(u_i)+V_3B_3^3(u_i)\big)\\
&=&c_{1,1}\alpha_1+c_{1,2}\alpha_2+\sum_{i = 1}^nA_{i,1}^\top\big(V_0B_0^3(u_i)+V_0B_1^3(u_i)+V_3B_2^3(u_i)+V_3B_3^3(u_i)\big)\tag{12}
\end{eqnarray}

式 $(12)$ で $c_{j,k}=\displaystyle\sum_{i = 1}^nA_{i,j}^\top A_{i,k}$ とおきました。
式 $(12)$ を式 $(11)$ に代入します。

\begin{eqnarray}
&&c_{1,1}\alpha_1+c_{1,2}\alpha_2+\sum_{i = 1}^nA_{i,1}^\top\big(V_0B_0^3(u_i)+V_0B_1^3(u_i)+V_3B_2^3(u_i)+V_3B_3^3(u_i)\big)=\sum_{i = 1}^nA_{i,1}^\top d_i\\
&&\Rightarrow c_{1,1}\alpha_1+c_{1,2}\alpha_2=\underbrace{\sum_{i = 1}^nA_{i,1}^\top\left(d_i - \left(V_0B_0^3(u_i)+V_0B_1^3(u_i)+V_3B_2^3(u_i)+V_3B_3^3(u_i)\right)\right)}_{=:X_1}\\
&&\Rightarrow c_{1,1}\alpha_1+c_{1,2}\alpha_2=X_1\tag{13}
\end{eqnarray}

式 $(13)$ で $X_j=\displaystyle\sum_{i = 1}^nA_{i,j}^\top\left(d_i - \left(V_0B_0^3(u_i)+V_0B_1^3(u_i)+V_3B_2^3(u_i)+V_3B_3^3(u_i)\right)\right)$ とおきました。
同様に式 $(8)$ は以下のようになります。

\begin{eqnarray}
c_{2,1}\alpha_1+c_{2,2}\alpha_2=X_2\tag{14}
\end{eqnarray}

式 $(13),(14)$ を変形していきます。

\begin{eqnarray}
&&\left\{
\begin{array}{l}
c_{1,1}\alpha_1+c_{1,2}\alpha_2=X_1\\
c_{2,1}\alpha_1+c_{2,2}\alpha_2=X_2
\end{array}
\right.\\
&&\Rightarrow
\begin{pmatrix}c_{1,1} & c_{1,2}\\ c_{2,1} & c_{2,2}\end{pmatrix}
\begin{pmatrix}\alpha_1\\ \alpha_2\end{pmatrix}
=
\begin{pmatrix}X_1\\ X_2\end{pmatrix}\\
&&\Rightarrow
\begin{pmatrix}\alpha_1\\ \alpha_2\end{pmatrix}
=
\begin{pmatrix}c_{1,1} & c_{1,2}\\ c_{2,1} & c_{2,2}\end{pmatrix}^{-1}
\begin{pmatrix}X_1\\ X_2\end{pmatrix}\\
&&\Rightarrow
\begin{pmatrix}\alpha_1\\ \alpha_2\end{pmatrix}
=
\dfrac{1}{\det\begin{pmatrix}c_{1,1} & c_{1,2}\\ c_{2,1} & c_{2,2}\end{pmatrix}}
\begin{pmatrix}c_{2,2} & -c_{1,2}\\ -c_{2,1} & c_{1,1}\end{pmatrix}
\begin{pmatrix}X_1\\ X_2\end{pmatrix}\\
&&\Rightarrow
\begin{pmatrix}\alpha_1\\ \alpha_2\end{pmatrix}
=
\dfrac{1}{\det\begin{pmatrix}c_{1,1} & c_{1,2}\\ c_{2,1} & c_{2,2}\end{pmatrix}}
\begin{pmatrix}c_{2,2} X_1 - c_{1,2}X_2\\ -c_{2,1}X_1 + c_{1,1}X_2\end{pmatrix}\\
&&\Rightarrow
\begin{pmatrix}\alpha_1\\ \alpha_2\end{pmatrix}
=
\dfrac{1}{\det\begin{pmatrix}c_{1,1} & c_{1,2}\\ c_{2,1} & c_{2,2}\end{pmatrix}}
\begin{pmatrix}
\det \begin{pmatrix} X_1 & c_{1,2}\\ X_2 & c_{2,2}\end{pmatrix}\\
\det \begin{pmatrix} c_{1,1} & X_1\\ c_{2,1} & X_2\end{pmatrix}
\end{pmatrix}\\
&&\Rightarrow
\begin{pmatrix}\alpha_1\\ \alpha_2\end{pmatrix}
=
\dfrac{1}{\det\begin{pmatrix}C_1 & C_2\end{pmatrix}}
\begin{pmatrix}
\det \begin{pmatrix} {\mathscr X} & C_2\end{pmatrix}\\
\det \begin{pmatrix} C_1 & {\mathscr X}\end{pmatrix}
\end{pmatrix}\\
&&\Rightarrow\left\{
\begin{array}{l}
\alpha_1=\dfrac{\det \begin{pmatrix} {\mathscr X} & C_2\end{pmatrix}}{\det\begin{pmatrix}C_1 & C_2\end{pmatrix}}\\
\alpha_2=\dfrac{\det \begin{pmatrix} C_1 & {\mathscr X}\end{pmatrix}}{\det\begin{pmatrix}C_1 & C_2\end{pmatrix}}\tag{15}
\end{array}
\right.
\end{eqnarray}

式 $(15)$ で $C_1=\begin{pmatrix}c_{1,1}\ c_{2,1}\end{pmatrix},\ C_2=\begin{pmatrix}c_{1,2}\ c_{2,2}\end{pmatrix},\ {\mathscr X}=\begin{pmatrix}X_1\ X_2\end{pmatrix}$ とおきました。
式 $(15)$より、$\alpha_1,\alpha_2$ が求まったので、式 $(5)$ より、$V_1,V_2$ が定まります。
$V_0,V_1,V_2,V_3$ が定まったので、3次ベジェ曲線が一旦求まったことになります。

⑤：求めたベジェ曲線と点列との最大距離を求め、指定誤差と比較する。

3次ベジェ曲線と点列との最大距離(以下、最大誤差と書く)を求めます。(ソースコードによると距離の2乗を求めていますが、元論文に合わせます。)

\newcommand{\argmax}{\mathop{\rm arg~max}\limits}\begin{eqnarray}
\max_i\ ||d_i-Q(u_i)||\hspace{30px}(i=2,\ldots n-1)\tag{16}
\end{eqnarray}

式 $(16)$ の値である最大誤差が指定した誤差 $\epsilon$ より小さければ、合格です。
点群の端点は3次ベジェ曲線と位置ベクトルが一致するので、調べる必要はありません。
そのため、式 $(9)$ の $i$ は端点を除いた $2$ から $n-1$ を対象としています。

⑥：粗めの誤差より最大誤差が小さい場合

ソースコードによると粗めの誤差 $\Psi=2\epsilon$ であり、最大誤差が $\Phi$ より小さい場合、本処理を行います。
(ソースコードによると粗めの誤差は指定誤差の4倍ですがそれは距離の2章で比較しているためであり、元論文では距離で比較しているので4倍の平方根である2倍にします。)
なぜ粗めの誤差で比較するのかというと、点列に問題があるのではなく、点列に対するパラメータの決め方が荒いため最大誤差が指定誤差の範囲おさまらなかったのだと考えます。

ある点 $P\in{\mathbb R}^2$ から曲線へ下した垂線と曲線の交点が $Q(t) であるとき、以下の式が成り立ちます。

\begin{eqnarray}
(Q(t)-P)^\top Q'(t)=0\tag{17}
\end{eqnarray}

式 $(17)$ を満たすような $t$ を求めるのですが、解析に解くのは困難なのでニュートン法で解きます。
式を簡単にするため、$f(t)=(Q(t)-P)^\top Q'(t),\ Q_1(t)=Q(t)-P,\ Q_2=Q'(t)$ とおくと、式 $(17)$ は以下のように書けます。

\begin{eqnarray}
f(t)=Q_1(t)^\top Q_2(t)=0\tag{18}
\end{eqnarray}

式 [tex:(18)] をニュートン法で解くには、

\begin{eqnarray}
t\leftarrow t - \dfrac{f(t)}{f'(t)}\tag{19}
\end{eqnarray}

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

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

\begin{eqnarray}
\dfrac{{\rm d}f}{{\rm d}t}&=&\dfrac{{\rm d}}{{\rm d}t} (Q(t)-P)^\top Q'(t)\\
&=& \dfrac{{\rm d}}{{\rm d}t}\left(Q(t)^\top Q'(t)-P^\top Q'(t)\right)\\
&=& \dfrac{{\rm d}}{{\rm d}t}Q(t)^\top Q'(t)-\dfrac{{\rm d}}{{\rm d}t}P^\top Q'(t)\\
&=& Q(t)'^\top Q'(t)+Q(t)^\top Q''(t)-P^\top Q''(t)\\
&=& || Q'(t) ||^2+(Q(t)-P)^\top Q''(t)\tag{20}
\end{eqnarray}

式 $(20)$ を式 $(19)$ に代入します。

\begin{eqnarray}
t\leftarrow t - \dfrac{(Q(t)-P)^\top Q'(t)}{ || Q'(t) ||^2+(Q(t)-P)^\top Q''(t)}\tag{21}
\end{eqnarray}

全ての $d_i$ について式 $(21)$ を一度だけ用いて、$u_i$ を更新します。(式 $(21)$ の $t$ を $u_i$ に置き換えてください。)
そして、再度 $\alpha_1,\alpha_2$ を求め、最大誤差を求めます。この最大誤差が $\epsilon$ 未満になれば、合格です。
ソースコードによると最大4回この処理 ($\{u_i\}_{i=1}^n$ と $\alpha_1,\alpha_2$ の更新) を行います。

特定の $\alpha_1,\alpha_2$ に対して、ニュートン法を一度だけ実行して $\{u_i\}_{i=1}^n$ を更新している理由ですが、
1回目のニュートン法では、参照する②で求めた $\alpha_1,\alpha_2$ は粗いといっても求めるべき値に近く、
ニュートン法は解に近いと少ない繰り返し回数で解に収束することが知られているため1回という少ない回数にしているものだと思われます。
2～4回目のニュートン法では、$\alpha_1,\alpha_2$ はさらに求めるべき値に近づいているので、1回で十分だと思われます。

⑦：粗めの誤差より最大誤差が大きい場合

最大誤差が $\Psi$ 以上の場合または最大誤差が $\Psi$ 未満であったが、最大誤差が $\epsilon$ 未満にならなかった場合の処理です。
最大誤差の点で点群を2分割します。
その際、分割点における単位接ベクトルを求めます。分割点の添え字を $i$ とすると、単位接ベクトルは次のように求まります。

\begin{eqnarray}
\hat{t}=\dfrac{d_{i - 1} - d_{i + 1}}{||d_{i - 1} - d_{i + 1}||}\tag{22}
\end{eqnarray}

以下のように分割した点群と端点の単位接ベクトルを用いて、3次ベジェ曲線をフィットさせます。
$1$ から $i$ までの点列と始点の単位接ベクトル $\hat{t}_1$ と終点の単位接ベクトル $\hat{t}$ を用いて、②に戻り処理を行います。
$i$ から $n$ までの点列と始点の単位接ベクトル $-\hat{t}$ と終点の単位接ベクトル $\hat{t}_2$ を用いて、②に戻り処理を行います。
分割点において向きだけ異なる単位接ベクトルを利用していますが、
これは分割点において2つの3次ベジェ曲線をG1連続にするためです。

これで終了です。お疲れさまでした。

最後に

本アルゴリズムは3次ベジェ曲線を接続するとき、(筆者が考えるには)G1連続で十分であるので、
先に接ベクトルを求めて、点列に対応するパラメータと制御点を更新して、フィットする3次ベジェ曲線を見つけるというものでした。
非常によくできたアルゴリズムだなと感心しました。

本アルゴリズムに対応した FitCurves.c のJavaScript への移植はこちらに書いております。

C0連続のみにしたい場合は、点群をその特徴点で分割して本アルゴリズムを適用すればよいと思います。
誤差を点群ごとに決めても面白いなと思っています。
点群のある点における接ベクトルの決定方法については、アイデアがあるので良い実験結果が得られれば、また記事を書きたいと考えております。

元論文の気づいた誤植を書いておきます。
p620 下 - $\alpha_1,\alpha_2$ になっていない。
p621 下 - $Q_1,Q_2$ の説明が抜けている。