補間: B-スプラインと平滑化フィッティング

概要

B-スプライン (Basis spline) は、低次の多項式を節点 (knot) でなめらかに接続した区分多項式を、 局所サポートを持つ基底関数 $N_{i,p}(t)$ の線形結合として表す枠組み。 曲線は制御点 $P_i$ の重み付き和 $C(t) = \sum_i N_{i,p}(t)\, P_i$ で書ける。

多項式補間に対する主な利点は次の 3 つ:

• 局所性: 各基底 $N_{i,p}$ は $p+1$ 個の節点区間だけで非零。 制御点を 1 つ動かしても曲線の近傍だけが変わり、全体には波及しない。
• 凸包性: 基底は非負で総和が $1$ (1 の分割) なので、曲線は各点で制御点の凸包に収まる。 これにより形状が暴れにくく、数値的にも安定。
• 任意次数: 次数 $p$ を選べる。等間隔データで高次多項式補間が起こす Runge 現象 (端での激しい振動) を、 低次の区分多項式で回避できる。

関連 API: bsplineCoefficients, bsplineEvaluate, bsplineRegression, splineFit。

Cox-de Boor 基底

次数 $p$ の B-スプライン基底関数は、節点ベクトル $\{t_0 \le t_1 \le \dots \le t_m\}$ 上でCox-de Boor の漸化式により再帰的に定義される。 次数 0 (区分定数) から出発する:

$$N_{i,0}(t) = \begin{cases} 1 & t_i \le t < t_{i+1} \\ 0 & \text{それ以外} \end{cases}$$

次数 $p \ge 1$ では、隣り合う 2 つの低次基底を $t$ の 1 次関数で重み付けして合成する:

$$N_{i,p}(t) = \frac{t - t_i}{t_{i+p} - t_i}\, N_{i,p-1}(t) + \frac{t_{i+p+1} - t}{t_{i+p+1} - t_{i+1}}\, N_{i+1,p-1}(t)$$

分母が $0$ になる項 (重複節点) は、その項全体を $0$ とみなす規約で扱う。 sangi の bsplineBasis はこの漸化式を直接実装し、分母が機械イプシロン以下の場合に当該項を落とす。

サポート区間

$N_{i,p}(t)$ は区間 $[t_i,\ t_{i+p+1})$ の外で恒等的に $0$ となる。 すなわち各基底は高々 $p+1$ 個の節点区間にしか影響しない。これが局所性の根拠であり、 評価時には非零の基底だけを考えればよいので計算量が抑えられる。

クランプ均等節点ベクトル

端点を曲線が確実に通るようにするため、両端の節点を多重度 $p+1$ で重ねた クランプ (clamped) 節点ベクトルを用いる。区間 $[a, b]$ 上で基底数 $n$ のとき、 全節点数は $n + p + 1$ で:

$$\underbrace{a, \dots, a}_{p+1},\ \underbrace{t_{p+1}, \dots, t_{n-1}}_{\text{内部節点}},\ \underbrace{b, \dots, b}_{p+1}$$

sangi の uniformKnots は内部節点を等間隔に置く。 補間 (bsplineCoefficients) では内部節点をデータ点の平均位置 (averaged knots) に取り、 コロケーション行列の良条件性を高める方式も用いる。

De Boor 評価

制御点 $P_i$ と節点ベクトルが与えられたとき、パラメータ $t$ での曲線値は

$$C(t) = \sum_{i} N_{i,p}(t)\, P_i$$

で定義される。これを基底関数を陽に計算せずに求めるのがDe Boor アルゴリズムである。 まず $t$ を含む節点区間 (knot span) $[t_k,\ t_{k+1})$ を特定し、 関係する $p+1$ 個の制御点 $P_{k-p}, \dots, P_k$ を初期値として三角形状に補間を繰り返す。

$r = 1, \dots, p$ の各段で、係数

$$\alpha = \frac{t - t_{k+1+j-p-1}}{t_{k+1+j-r} - t_{k+1+j-p-1}}$$

を用いて

$$d_j^{(r)} = (1 - \alpha)\, d_{j-1}^{(r-1)} + \alpha\, d_j^{(r-1)}$$

と更新し、最終段の $d_p^{(p)}$ が曲線値 $C(t)$ になる。

数値安定性

各更新は係数 $\alpha \in [0, 1]$ による凸結合 (convex combination) である。 値の補間が常に 2 点の内分で行われるため、係数が増幅されず丸め誤差が蓄積しにくい。 基底を再帰展開して総和を取る素朴な方法より、De Boor の三角アルゴリズムは数値的に頑健である。 sangi の deBoor は分母が機械イプシロン以下のとき $\alpha = 0$ とし、重複節点でも破綻しないようにしている。

B-スプライン補間と回帰

補間 (interpolation)

データ点 $(x_i, y_i)$ ($i = 0, \dots, n-1$) をすべて厳密に通す係数 $c_j$ を求める。 基底行列 (コロケーション行列) を $B_{ij} = N_{j,p}(x_i)$ とすると、補間条件 $\sum_j c_j N_{j,p}(x_i) = y_i$ は 正方の線形系

$$B\, c = y$$

になる。$B$ は基底の局所性から帯行列となり、効率よく解ける。 sangi の bsplineCoefficients はクランプ節点でこの系を LU 分解で解く。 さらに端点行 ($x_0$ と $x_{n-1}$) は $N_{0,p}(x_0) = 1$, $N_{n-1,p}(x_{n-1}) = 1$ という事実を使って明示的に設定し、 端点での数値誤差を防いでいる。

回帰 (regression)

基底数 $M$ をデータ点数 $n$ より少なく取ると、$B$ は $n \times M$ の縦長行列になり厳密解は一般に存在しない。 そこで残差平方和 $\|y - Bc\|^2$ を最小化する最小二乗解を、正規方程式

$$B^\top B\, c = B^\top y$$

から求める。$B^\top B$ は $M \times M$ の対称半正定値行列である。 基底数が少ないほど曲線は滑らかになり、データのノイズを平均化できる。

sangi の bsplineRegression は正規方程式を構成したのち、対角に小さな正則化項

$$(B^\top B + \lambda I)\, c = B^\top y, \qquad \lambda = 10^{-10}$$

を加えて Gauss 消去で解く。$\lambda$ は $B^\top B$ がほぼ特異 (基底が一部のデータ範囲で励起されない等) なときの数値的破綻を防ぐリッジ項であり、 フィットの形状にはほとんど影響しない大きさに取る。

平滑化フィッティングと平滑化パラメータ $s$

補間 ($s=0$) と強い平滑化の間を連続的に橋渡しするのが平滑化スプラインである。 重み $w_i$ 付きの残差平方和を、平滑化パラメータ $s$ を上限として制約する:

$$\sum_{i} w_i\, \bigl(y_i - S(x_i)\bigr)^2 \le s$$

• $s = 0$: 全データ点を厳密に通る補間。
• $s$ 大: 残差を許容して節点を減らし、滑らかな曲線。ノイズの影響を抑える。

$s$ は残差 (データへの忠実さ) と滑らかさのトレードオフを司る。 sangi の splineFit は、内部節点を少しずつ増やしながら各回で重み付き最小二乗フィットを行い、 残差が $s$ 以下になった時点で停止する (過適合を避けるため、残差が悪化に転じたら早期終了する)。 $s$ を負に指定すると、Dierckx の推奨に従って $s = n$ (データ点数) を既定値として自動選択する。

各候補節点数での最小二乗は、重み行列 $W = \mathrm{diag}(w_i)$ を用いた重み付き正規方程式

$$(B^\top W B)\, c = B^\top W y$$

として解く。実装は $W^{1/2}B$ と $W^{1/2}y$ を作ってから正規方程式を組み、LU 分解で求解する。

パラメトリック曲線・サーフェススプライン

弦長パラメトリック曲線

$y$ が $x$ の関数でない閉曲線や折れ曲がった軌跡は、$y = f(x)$ の形では表せない。 この場合、共通のパラメータ $t$ に対して $x(t)$ と $y(t)$ を別々のスプラインでフィットする。

パラメータには弦長 (chord length) パラメータ化を用いる。 隣接データ点間のユークリッド距離を累積し、全長で正規化する:

$$t_0 = 0, \qquad t_i = t_{i-1} + \sqrt{(x_i - x_{i-1})^2 + (y_i - y_{i-1})^2}, \qquad t_i \leftarrow t_i / t_{n-1}$$

これにより点が密な区間ではパラメータがゆっくり進み、曲線の形状によらず偏りの少ないフィットが得られる。 sangi の parametricSplineFit は $x(t)$, $y(t)$ をそれぞれ splineFit でフィットし、 評価は両者の値の組 $(x(t), y(t))$ として返す。

2D テンソル積サーフェススプライン

格子状データ $z_{ij} = f(x_i, y_j)$ をなめらかな曲面で表すには、 2 方向の B-スプライン基底の積 (テンソル積) を用いる:

$$S(x, y) = \sum_i \sum_j c_{ij}\, N_i(x)\, N_j(y)$$

係数 $c_{ij}$ は2 段階の 1 次元求解で得られる。 まず各 $x$ 行について $y$ 方向の B-スプライン係数を求め (中間係数 $\alpha_{ij}$)、 次にその各 $y$ 列について $x$ 方向の B-スプライン係数を求める。 $x$ 方向・$y$ 方向それぞれにクランプ節点ベクトルを作り、各 1 次元系を LU 分解で解く。 sangi の surfaceSplineFit がこの分離求解を実装し、評価 surfaceSplineEval は 非零の $x$ 方向基底だけをループして 2 重和を計算する。

画像補間

画像のリサンプリング (拡大・縮小・回転) は、整数格子上の画素値 $g[r][c]$ を 非整数の評価点 $(x, y)$ で再構成する問題である。 再構成は補間カーネル $K$ による畳み込み $\hat g(x, y) = \sum_{m, n} K(x - n)\, K(y - m)\, g[m][n]$ で表され、 カーネルの台 (サポート) の広さと形が再構成品質を決める。 sangi は次の 3 種を提供する。

最近傍 (nearest neighbor)

評価点をもっとも近い格子点に丸め、その画素値をそのまま返す (nearestNeighbor2D)。 台は 1 画素分。最速だが段差状の不連続が出るため、ピクセルアートの整数倍拡大など、にじみを避けたい用途に向く。

双 3 次 (bicubic, Catmull-Rom カーネル)

$4 \times 4 = 16$ 個の近傍画素を 3 次多項式カーネルで重み付けする (bicubicInterpolate)。 sangi はCatmull-Rom カーネル ($a = -1/2$) を採用する:

$$K(t) = \begin{cases} \tfrac{3}{2}|t|^3 - \tfrac{5}{2}|t|^2 + 1 & |t| \le 1 \\[4pt] -\tfrac{1}{2}|t|^3 + \tfrac{5}{2}|t|^2 - 4|t| + 2 & 1 < |t| < 2 \\[4pt] 0 & |t| \ge 2 \end{cases}$$

このカーネルは台が $[-2, 2]$ で、格子点で補間条件 $K(0)=1$, $K(\pm 1)=K(\pm 2)=0$ を満たす。 なめらかでエッジも比較的保たれるため、一般的なリサンプリングの既定として広く使われる。

Lanczos リサンプリング

$\mathrm{sinc}$ 関数を窓掛けしたLanczos カーネルで、もっとも鮮鋭な再構成が得られる (lanczosInterpolate2D):

$$K(t) = \begin{cases} \mathrm{sinc}(t)\, \mathrm{sinc}(t/a) & |t| < a \\ 0 & |t| \ge a \end{cases}, \qquad \mathrm{sinc}(t) = \frac{\sin(\pi t)}{\pi t}$$

$a$ は窓幅で、$a = 2$ なら $4 \times 4$ 窓、$a = 3$ (既定) なら $6 \times 6$ 窓を使う。 台が広いほど理想的なローパスフィルタに近づき高品質だが、計算量が増え、輪郭近傍にリンギングが出ることがある。 sangi の実装は境界付近のエネルギー保存のために重みの総和で正規化する。

カーネルの台と品質の関係は、台が広いほど周波数特性が理想に近づく一方で計算コストとリンギングが増す、 という一貫したトレードオフとして整理できる。

比較表

補間 vs 回帰 vs 平滑化

計算量

大まかな指針: データに誤差がなければ補間、ノイズがあれば回帰または平滑化 ($s$ を調整)、 閉曲線や軌跡なら弦長パラメトリックフィット、格子曲面ならテンソル積、 画像リサンプリングは品質と速度で最近傍 / 双 3 次 / Lanczos を選ぶ。

参考文献

• de Boor, C. (2001). A Practical Guide to Splines. Revised ed. Springer.
• Dierckx, P. (1993). Curve and Surface Fitting with Splines. Oxford University Press.
• Cox, M. G. (1972). "The numerical evaluation of B-splines". Journal of the Institute of Mathematics and its Applications, 10(2), 134–149.
• Piegl, L., & Tiller, W. (1997). The NURBS Book. 2nd ed. Springer.
• Catmull, E., & Rom, R. (1974). "A class of local interpolating splines". In Computer Aided Geometric Design (pp. 317–326). Academic Press.
• Duchon, C. E. (1979). "Lanczos filtering in one and two dimensions". Journal of Applied Meteorology, 18(8), 1016–1022.