軸方向に単調なBezier曲線の最小二乗法

最小二乗法を用いて, 指定点を通るようなBezier曲線を求めたい. はじめに, Bezier曲線のパラメータの関係で, 例えばx, y平面を考える場合は, xまたはyについて単調増加(あるいは単調減少)である必要がある. 最後にPythonによって実装を書いておいた.

間違っている可能性があるので注意.

[目次]

Bezier曲線の定義
Bezier曲線による最小二乗法
Pythonで実装

Bezier曲線の定義

空間上に点 $\boldsymbol{P}_0, \boldsymbol{P}_n$ があるとする. 制御点 $\boldsymbol{P}_1, \cdots, \boldsymbol{P}_{n-1}$ とパラメータ $t (0 \leq t \leq 1)$ によって, $n$ 次元Bezier曲線は次のように定義される;

$\displaystyle \boldsymbol{R} \left(t \right) = \sum^{n}_{i=0} B^{n}_{i} \left(t \right) \boldsymbol{P}_{i} \tag{1}$

ここに $B^{n}_{i} \left(t \right)$ は次式で定義されるBernstein関数である;

$\displaystyle B^{n}_{i} \left(t \right) = \frac{n!}{i! \left(n-i \right)!} t^i \left(1-t \right)^{n-i} \tag{2}$

以降は点 $\boldsymbol{P}_0, \boldsymbol{P}_n$ と制御点 $\boldsymbol{P}_1, \cdots, \boldsymbol{P}_{n-1}$ とをまとめて制御点と呼ぶことにする( $\boldsymbol{R} \left(t \right)$ の計算から, その方が便利なことがわかる). なお $t=0, t=1$ を $\boldsymbol{R} \left(t \right)$ に代入してみれば明らかなように, Bezier曲線は必ず点 $\boldsymbol{P}_0, \boldsymbol{P}_n$ を通る.

Bezier曲線による最小二乗法

空間上にBezier曲線が通ることを期待したい点 $\boldsymbol{Q}_j = \left(x_j, y_j, z_j \right)^{T} \left(j \in \left\{0, 1, \cdots, m \right\}\right)$ があり, 点 $\boldsymbol{Q}_j$ 各々に対してBezier曲線のパラメータ $t_j$ がわかっているとする. このとき, 点 $\boldsymbol{Q}_j$ を通るような $n$ 次元Bezier曲線の制御点 $\boldsymbol{P}_i = \left(x_i, y_i, z_i \right)^{T} \left(i \in \left\{0, 1, \cdots, n \right\}\right)$ を見つけたい. ただし上側添字 $T$ は転置を意味する. 最小二乗法を考えれば, 点 $\boldsymbol{Q}_j$ について

$\displaystyle E \left( \boldsymbol{Q}_j \right) = \left(\boldsymbol{Q}_j - \boldsymbol{R} \left(t_j, \boldsymbol{P} \right) \right) \cdot \left(\boldsymbol{Q}_j - \boldsymbol{R} \left(t_j, \boldsymbol{P} \right) \right) \tag{3}$

を最小にするような制御点 $\boldsymbol{P}_i$ を見つけることになる. ただし, $\boldsymbol{P}=\left\{ \boldsymbol{P}_0, \boldsymbol{P}_1, \cdots, \boldsymbol{P}_n \right\}$ である. 今の問題では制御点 $\boldsymbol{P}=\left\{ \boldsymbol{P}_0, \boldsymbol{P}_1, \cdots, \boldsymbol{P}_n \right\}$ がパラメータのようなものなので, (1)で与えられるBezier曲線を $\boldsymbol{P}=\left\{ \boldsymbol{P}_0, \boldsymbol{P}_1, \cdots, \boldsymbol{P}_n \right\}$ の関数とした.

ここで問題は, 仮定から明らかなように, 点 $\boldsymbol{Q}_j$ とパラメータ $t_j$ があらかじめわかっていなければならないことである. 点 $\boldsymbol{Q}_j$ だけが与えられ, 点 $\boldsymbol{Q}_j$ を通るBezier曲線を求めたい場合, $t_j$ がわかっている場合は少ない(と筆者は思う)ので, この方法は(筆者にとっては)微妙である......( $t=0,t=1$ の場合についてはわかっているがそれ以外がわからない.)

そこで, 点 $\boldsymbol{Q}_j$ のパラメータ $x, y, z$ のいずれか1つの単調性を仮定する. 単調変化するパラメータの座標値を0から1の範囲に規格化することでパラメータを代用することができる. 例えばパラメータxが単調変化する場合は, 規格化することで, $y = R\left(x\right)$ となる！

以下では簡単のために, $x, y$ 平面で点 $\boldsymbol{Q}_j$ のパラメータ $x$ について単調性を仮定する. パラメータ $x$ をBezier曲線のパラメータ $t$ とみなすために, 先に点 $\boldsymbol{Q}_j$ の $x$ 座標の最大値 $x_{max}$ および最小値 $x_{min}$ を用いて, Min-Max Normalizationを行っておく;

$\displaystyle x_{j} \longleftarrow \frac{x_{j} - x_{min}}{x_{max} - x_{min}} \tag{4}$

ただしこの方法では $\boldsymbol{Q}_0$ が $t_0$ 、 $\boldsymbol{Q}_m$ が $t_m$ に対応する*1. そこで座標系の最大値と最小値を用いても良い； $x$ 座標が $[x_{min}, x_{max}]$ の範囲だと思って式(4)を使えば良い.

これによってパラメータを代用できる(Bezier曲線を求めたら, 上の計算を $x_j$ について計算すれば, 元の座標系に戻すことができる). また1つの軸をパラメータに見立てることで, 制御点 $\boldsymbol{P}_i$ はベクトルで計算する必要がなくなり, 制御値 $P_i$ ( $y$ 座標)を求めるだけでよくなるため, 式(3)は次のようになる;

$\displaystyle E \left( \boldsymbol{Q}_j \right) = \left(y_j - R \left(x_j, \boldsymbol{P} \right) \right)^2 \tag{5}$
ただし, $\boldsymbol{P}=\left\{ P_0, P_1, \cdots, P_n \right\}$ である. 式(5)を制御値 $P_k \left(k \in \left\{0, 1, \cdots, n \right\}\right)$ で微分する. $R \left(x_j, \boldsymbol{P} \right)$ はスカラー $P_k$ の一次関数なので瞬時に計算できて
$\displaystyle \frac{\partial E \left( \boldsymbol{Q}_j \right) }{\partial P_k} = - 2 B^{n}_{k} \left(x_j \right) \left(y_j - R \left(x_j, \boldsymbol{P} \right) \right) \tag{6}$

$j$ について足し合わせる. これが0となるような $\boldsymbol{P}$ を探していることを思い出す;
$\begin{align*} \frac{1}{2} \sum^{m}_{j=0} \frac{\partial E \left( \boldsymbol{Q}_j \right) }{\partial P_k} &= - \sum^{m}_{j=0} B^{n}_{k} \left(x_j \right) \left(y_j - R \left(x_j , \boldsymbol{P} \right) \right) \\\\ &= - \sum^{m}_{j=0} B^{n}_{k} \left(x_j \right) \left(y_j - \sum^{n}_{i=0} B^{n}_{i} \left(x_j \right) P_{i} \right) \tag{7} \\\\ &= 0 \end{align*}$

展開して制御値 $P_i$ について整理すれば,

$\displaystyle P_{0} \sum^{m}_{j=0} B^{n}_{0} \left( x_j \right) B^{n}_{k} \left( x_j \right) + P_{1} \sum^{m}_{j=0} B^{n}_{1} \left( x_j \right) B^{n}_{k} \left( x_j \right) + \ldots \\\\ \displaystyle + P_{n} \sum^{m}_{j=0} B^{n}_{n} \left( x_j \right) B^{n}_{k} \left( x_j \right) = \sum^{m}_{j=0} B^{n}_{k} \left(x_j \right) y_j \tag{8}$

頑張って整理すると次のようになる;

$\begin{align*} \left( \begin{array}{cccc} B^{n}_{0} \left( x_0 \right) & B^{n}_{1} \left( x_0 \right) & \cdots & B^{n}_{n} \left( x_0 \right) \\ B^{n}_{0} \left( x_1 \right) & B^{n}_{1} \left( x_1 \right) & \cdots & B^{n}_{n} \left( x_1 \right) \\ \vdots & \vdots & \ddots & \vdots \\ B^{n}_{0} \left( x_m \right) & B^{n}_{1} \left( x_m \right) & \cdots & B^{n}_{n} \left( x_m \right) \end{array} \right) \left( \begin{array}{cccc} P_{0} \\ P_{1} \\ \vdots \\ P_{n} \end{array} \right) = \left( \begin{array}{c} y_0 \\ y_1 \\ \vdots \\ y_m \end{array} \right) \tag{10} \end{align*}$

この方程式を解けば, 制御値 $\boldsymbol{P} = \left(P_0, P_1, \cdots, P_n \right)^{T}$ を求めることができる(逆行列を作用させるだけで良い).

Pythonで実装

これをPythonで実装した. 環境は下の通り.

Python==3.8.3
numpy==1.19.1
scipy==1.5.2
matplotlib==3.3.1

from scipy.special import comb
import numpy as np
import matplotlib.pyplot as plt

def bernstein(n, i, t):
    return comb(n, i) * t**i * (1 - t)**(n-i)

def bezier(n, t, P):
    return np.dot([bernstein(n, i, t) for i in range(n + 1)], P)

if __name__ == '__main__':
    # 通りたい点Q
    Q = np.array([(0, 0), (1, 2), (3, 0), (4, 4)])
    x = Q[:, 0]
    y = Q[:, 1]

    # 今回はQのx座標が単調増加しているので,
    # xをパラメータtに変換する.
    x_max = x.max()
    x_min = x.min()
    t = (x - x_min) / (x_max - x_min)

    # Bezier曲線の次元.
    n = 3

    # 式(10)の方程式を解く.
    # Qによっては逆行列が存在しないので,
    # np.linalg.solve(方程式)やnp.linalg.inv(逆行列)では計算できない.
    # そこでnp.linalg.pinv(模擬逆行列)で強引に計算している.
    A = np.array([[bernstein(n, i, tj) for i in range(n+1)] for tj in t])
    P = np.dot(np.linalg.pinv(A), y)  # 制御値

    # 求めたBezier曲線
    x_lsm = np.linspace(0, 1, 100) * (x_max - x_min) + x_min
    y_lsm = np.array([bezier(n, t, P) for t in np.linspace(0, 1, 100)])

    plt.figure()
    plt.plot(x_lsm, y_lsm, label="estimated")  # Bezier曲線
    plt.plot(x, y, '--x', label="Q")  # 入力点
    plt.plot(x, P, '--o', label="P")  # 制御値
    plt.legend()
    plt.show()