N次ベジェ曲線によるカーブフィッティング

はじめに

先日，ベジェ曲線上の点ををなるべく短いコード計算するだけの記事を投稿した．
【Python】なるべく短くベジェ曲線

今回は，任意の $m$ 個のデータ点 $\boldsymbol{p}_i=(x_i, y_i)$ (もっと多次元でも可)に１本の $N$ 次ベジェ曲線をフィッティングしてみようと思う．
既存の記事は，3次ベジェ曲線に限定していたり，端点を固定していたり，そもそも説明と実装が対応しているかも怪しかったりするので．

ベジェ曲線はデザインに多用されることからも分かるように非常に表現力が高く，癖も少ないので，多項式フィッティングのように「外挿しようとしたらフィッティング範囲外では変な値に発散していた...」「ぱっと見良さそうだけど誤差を計算すると波打ってる」みたいな現象は比較的少ないのではないかと思う（ただし外挿する場合は端点の $x$ 座標を固定してフィッティングしないと定義域外は使用できないので注意）．
欠点もないわけではなく，特に今回は1本のベジェ曲線でフィッティングするので直線の表現が苦手だったり，媒介変数で表現されるために任意の $x$ に対応する $y$ を知りたい場合は線形補間など何らかの手段で対応関係を特定しないといけなかったりする．
また，フィッティング手法によると思われる問題として，外れ値に妙に強く影響される点がある．そのため実データで解析を行う場合は事前にクレンジングが必要である．

フィッティングはだいたい最小二乗法の手法に倣って行う．

最小二乗法

$\boldsymbol{p_{i}}=(x_i, y_i)$ を $y_i\sim f(x_i)$ のようにフィッティングする場合，

e = \sum_i \left(y_i - f(x_i)\right)^2

が最小になるようにする．例えば $f(x)=ax+b$ とすると，

e = \sum_i \left(y_i - a x_i-b\right)^2 = \sum_i \left[a^2 x_i^2 - 2(y_i - b)ax_i + b^2 - 2y_ib + y_i^2 \right]

これが最小となるパラメータ $a, b$ を知りたい．それには上の式をパラメータで微分した値が $0$ になる $a, b$ を求めれば良いので，

\frac{\partial e}{\partial a} = \sum_i \left[2a x_i^2 - 2(y_i - b)x_i \right] = 0\\
\frac{\partial e}{\partial b} = \sum_i \left[ 2ax_i + 2b - 2y_i \right] = 0

$y_i$ を含む項を右辺に移動して，行列表示にすると，

\left(
\begin{array}\\
 \sum_i x_i^2 & \sum_i x_i\\
 \sum_i x_i & \sum_i 1
\end{array}
\right)
\left(
\begin{array}\\
 a\\
 b
\end{array}
\right)
=
\left(
\begin{array}\\
 \sum_i x_i y_i \\
 \sum_i y_i
\end{array}
\right)

であるから，

\left(
\begin{array}\\
 a\\
 b
\end{array}
\right)
=
\left(
\begin{array}\\
 \sum_i x_i^2 & \sum_i x_i\\
 \sum_i x_i & \sum_i 1
\end{array}
\right)^{-1}
\left(
\begin{array}\\
 \sum_i x_i y_i \\
 \sum_i y_i
\end{array}
\right)

を計算すればよい．

ベジェ曲線に最小二乗法（様のもの）を適用する

$N$ 次ベジェ曲線の定義は以下である．

\boldsymbol{P}(t)=\sum_{i=0}^{N}\boldsymbol{B}_i\left(\begin{array}{c}n\\ i\\ \end{array}\right)t^i(1-t)^{n-i}\\
\mathrm{where} \quad t \in [0, 1]

$\boldsymbol{B_i}=(X_i, Y_i)$ は制御点の座標である．
ベジェ曲線の定義は媒介変数 $t$ によって表現されるため，$x, y$ が出現しない．困った．

とりあえずこのまま計算してみよう．知りたいパラメータは $N$ 個の制御点の座標 $\boldsymbol{B_i}=(X_i, Y_i)$ である．

\left(\begin{array}{c}n\\ i\\ \end{array}\right)t^i(1-t)^{n-i} \quad(\equiv f(i, t))

は $\boldsymbol{P_i}$に依存しないので$\boldsymbol{F}(t)=(f(0, t),...,f(N-1, t))^T$と置き直して，

\begin{align}
e &= \sum_i \left(\boldsymbol{p}_i - \boldsymbol{P}(t)\right)^2\\
&= \sum_i \left(\boldsymbol{p}_i - \boldsymbol{B}\boldsymbol{F} \right)^2\\
&= \sum_i \left[\boldsymbol{p}_i^2 + (\boldsymbol{B}\boldsymbol{F})^2 - 2\boldsymbol{p}_i\boldsymbol{B}\boldsymbol{F} \right]
\end{align}

ここで$(\boldsymbol{B}\boldsymbol{F})^2$の$\boldsymbol{B}_i$に関係する部分を$(\boldsymbol{B}\boldsymbol{F})^2_i$と表現すると，

\begin{align}
(\boldsymbol{B}\boldsymbol{F})^2_i = \boldsymbol{B}_i^2f(i, t)^2 + 2\boldsymbol{B}_if(i, t)\sum_{j\neq i}\boldsymbol{B}_jf(j, t)
\end{align}

である．パラメータで微分して，

\begin{align}
\frac{\partial e}{\partial \boldsymbol{B}_j}
&= \sum_i \left[\frac{\partial}{\partial \boldsymbol{B_j}}(\boldsymbol{B}\boldsymbol{F})^2 - 2\boldsymbol{p}_if(j, t_i) \right]\\
&= \sum_i \left[2\boldsymbol{B}_jf(j, t_i)^2 + 2f(j, t_i)\sum_{k\neq j}\boldsymbol{B}_kf(k, t_i) - 2\boldsymbol{p}_if(j, t_i) \right]\\
&=0
\end{align}

$\boldsymbol{p}_i$ を含む項を右辺に移動して，行列表示にすると，

\left(
\begin{array}\\
 \sum_i f(0, t_i)^2 & \sum_i f(0, t_i)f(1, t_i) & \cdots & \sum_i f(0, t_i)f(N-1, t_i) \\
 \sum_i f(0, t_i)f(1, t_i) & \sum_i f(1, t_i)^2 & \cdots & \sum_i f(1, t_i)f(N-1, t_i) \\
\vdots &\vdots &\ddots &\vdots\\
 \sum_i f(N-1, t_i)f(0, t_i) & \sum_i f(N-1, t_i)f(1, t_i) & \cdots & \sum_i f(N-1, t_i)^2\\
\end{array}
\right)
\left(
\begin{array}\\
 \boldsymbol{B}_0\\
 \boldsymbol{B}_1\\
 \vdots \\
 \boldsymbol{B}_{N-1}
\end{array}
\right)
=
\left(
\begin{array}\\
 \sum_i \boldsymbol{p}_if(0, t_i) \\
 \sum_i \boldsymbol{p}_if(1, t_i)\\
 \vdots\\
 \sum_i \boldsymbol{p}_if(N-1, t_i)\\
\end{array}
\right)

であるから，

\left(
\begin{array}\\
 \boldsymbol{B}_0\\
 \boldsymbol{B}_1\\
 \vdots \\
 \boldsymbol{B}_{N-1}
\end{array}
\right)
=
\left(
\begin{array}\\
 \sum_i f(0, t_i)^2 & \sum_i f(0, t_i)f(1, t_i) & \cdots & \sum_i f(0, t_i)f(N-1, t_i) \\
 \sum_i f(0, t_i)f(1, t_i) & \sum_i f(1, t_i)^2 & \cdots & \sum_i f(1, t_i)f(N-1, t_i) \\
\vdots &\vdots &\ddots &\vdots\\
 \sum_i f(N-1, t_i)f(0, t_i) & \sum_i f(N-1, t_i)f(1, t_i) & \cdots & \sum_i f(N-1, t_i)^2\\
\end{array}
\right)^{-1}
\left(
\begin{array}\\
 \sum_i \boldsymbol{p}_if(0, t_i) \\
 \sum_i \boldsymbol{p}_if(1, t_i)\\
 \vdots\\
 \sum_i \boldsymbol{p}_if(N-1, t_i)\\
\end{array}
\right)

を解けば良い．
...のだが，$i$ についての和 $\sum_i$ の $t=t_i$ はどう扱えば良いのかという問題が残っていた．

何が正解かわからないので，仮置きした $\boldsymbol{B}$ から計算した $\boldsymbol{P}(t)$ のうち最もユークリッド距離の近いインデックスの値で代用してみる．
つまり，

t_i = \mathrm{argmin}(||\boldsymbol{p}_i-\boldsymbol{P}(t)||)

である．単純だが案外上手くいく．実際に確かめてみよう．

（$x$成分のみでという選択肢もある．）

t_i = \mathrm{argmin}(|x_i-\boldsymbol{P}(t)_x|)

Pythonで計算してみる

ベジェ曲線の計算には前回の記事のコードを少し書き換えて使う．

import numpy as np
import matplotlib.pyplot as plt

def F(n, t):
    J = np.vectorize(lambda i, t: (n/np.r_[1:i+1]-1).prod() * t**i * (1-t)**(n-1-i), signature='(),(m)->(m)')
    return J(np.r_[:n], t)

def Bezier(B, t):
    *_, n = B.shape
    return B @ F(n, t)

パラメータの準備

# ベジェ曲線のサンプル点
t = np.linspace(0, 1, 1000)

# 適当な(x, y)
x = np.linspace(-3, 3, 40)
y = np.sin(1.5*x) + .2*np.random.randn(*x.shape)

# 次数+1，とりあえずスタンダードに3次ベジェ
n = 4

# Bを初期化，変な値だと変にフィッティングされるので(x,y)から均等にとってみる（順番がバラバラの場合ソートしておいたほうが良いかも）
j = np.linspace(0, len(x)-1, n, True, dtype=int)
B = np.array([[x[i] for i in j], [y[i] for i in j]])

計算

for i in range(5): # 一回で収束しないことが多いので何回かイテレーション
    f = F(n, t)
    idx = ((Bezier(B, t).T[...,None] - np.r_[[x, y]])**2).sum(1).argmin(0) # argmin(pi-P(t))
    inv = np.linalg.pinv(f[:, idx] @ f[:, idx].T) # (疑似)逆行列
    # 成分ごとに更新
    x_ = inv @ (x * f[:, idx]).sum(1, keepdims=True)
    y_ = inv @ (y * f[:, idx]).sum(1, keepdims=True)
    B = np.c_[x_, y_].T # Bを更新

# プロット
plt.plot(*Bezier(B, t), c='r', lw=1)
plt.plot(*B, marker='o', ls='--', lw=.5, c='k', alpha=.5)
plt.scatter(x, y, s=5)
plt.show()

上手くフィットしている．

3次元も試してみる

t = np.linspace(0, 1, 1000)
x = np.linspace(-3, 3, 200)
y = np.sin(1.5*x) + .2*np.random.randn(*x.shape)
z = np.cos(1.5*x) + .2*np.random.randn(*x.shape) # z を追加．位相を90度ずらしたので螺旋になるはず．

n = 7 # 6次ベジェ
j = np.linspace(0, len(x)-1, n, True, dtype=int)
B = np.array([[x[i] for i in j], [y[i] for i in j], [z[i] for i in j]])

# 計算
for i in range(5):
    f = F(n, t)
    idx = ((Bezier(B, t).T[...,None] - np.r_[[x, y, z]])**2).sum(1).argmin(0)
    inv = np.linalg.pinv(f[:, idx] @ f[:, idx].T)
    x_ = inv @ (x * f[:, idx]).sum(1, keepdims=True)
    y_ = inv @ (y * f[:, idx]).sum(1, keepdims=True)
    z_ = inv @ (z * f[:, idx]).sum(1, keepdims=True)
    B = np.c_[x_, y_, z_].T

fig, ax = plt.subplots(subplot_kw={'projection': '3d'}) # 3D
ax.plot(*Bezier(B, t), c='r', lw=1)
ax.plot(*B, marker='o', ls='--', lw=.5, c='k', alpha=.5)
ax.scatter(x, y, z, s=5)
plt.show()

バッチリ．

陰関数もいけるのでは

# ベジェ曲線のサンプル点
t = np.linspace(0, 1, 1000)

# 円
a = np.linspace(0, 2*np.pi, 400)
x = np.cos(a)
y = np.sin(a)

n = 6 # 5次ベジェ

j = np.linspace(0, len(x)-1, n, True, dtype=int)
B = np.array([[x[i] for i in j], [y[i] for i in j]])

# 計算
for i in range(10):
    f = F(n, t)
    idx = ((Bezier(B, t).T[...,None] - np.r_[[x, y]])**2).sum(1).argmin(0) # argmin(pi-P(t))
    inv = np.linalg.pinv(f[:, idx] @ f[:, idx].T)
    x_ = inv @ (x * f[:, idx]).sum(1, keepdims=True)
    y_ = inv @ (y * f[:, idx]).sum(1, keepdims=True)
    B = np.c_[x_, y_].T

# プロット
plt.plot(*Bezier(B, t), c='r', lw=1)
plt.plot(*B, marker='o', ls='--', lw=.5, c='k', alpha=.5)
plt.scatter(x, y, s=5)
plt.axis("equal")
plt.show()

いけた．

おわりに

修論発表のスライドが作りたくない一心で新たな記事が誕生した．結局LaTeX書いてるんだよなあ．