先に投稿した『バイキュービック補間の1パス処理』の内容を、実装技術者向け入門編の一部としてまとめ直しました。

はじめに

画像処理におけるバイキュービック補間は、各画素の輝度値から画素と画素の間を、3次関数によりxおよびy方向に補間するものです。画素と画素の間の値を求めることで、画像をずらす（シフトする）ことや、拡大・縮小することができます。

バイキュービック畳み込みアルゴリズム

画像処理のバイキュービック補間では、各画素について補間計算をするのは大変なので、あらかじめ畳み込みカーネルを求めておき、入力画像の各画素に畳み込む畳み込みアルゴリズムが用いられます。
畳み込みアルゴリズムにおける入力および補間関数やデータ列の定義は以下の通りです。

	呼称	表記
1	入力関数	$f(x)$
2	入力データ列	$f(x_j)$
3	補間関数	$g(x)$
4	補間データ列

1と3は連続関数です。
通常の画像処理では、1は観測できません。2を入力として4を求めます。
1をサンプリングしたものが2です。3を再サンプリングしたものが4です。
サンプリングの間隔は一定ですが、再サンプリングの間隔は任意です。
$x$の座標系は任意です。

畳み込みカーネルはデルタ関数（$s=0$だけ$1$、他は$0$）上の点を補間する以下の3次関数の合成（区分関数）です。

$1 \le s \le 2$ については、点$(1,0)$、$(2,0)$を通る$W_0(s)$
$0 \le s \le 1$ については、点$(0,1)$、$(1,0)$を通る$W_1(s)$
$-1 \le s \le 0$ については、点$(-1,0)$、$(0,1)$を通る$W_2(s)$
$-2 \le s \le -1$ については、点$(-2,0)$、$(-1,0)$を通る$W_3(s)$
上記以外については$0$

ただし、$s$は畳み込みカーネルの座標系で；
$\quad s=(x-x_j)/h$
とします。$x_j$は$x$以下の直近のサンプリング位置、$h$はサンプリング間隔です。

畳み込みカーネルが滑らか（1次微分が連続）になるよう、以下の条件を追加します。（$a$は定数）

$W_0(s)$の1次微分$W_0'(s)$は、点$(1,a)$、$(2,0)$を通る
$W_1(s)$の1次微分$W_1'(s)$は、点$(0,0)$、$(1,a)$を通る
$W_2(s)$の1次微分$W_2'(s)$は、点$(-1,-a)$、$(0,0)$を通る
$W_3(s)$の1次微分$W_3'(s)$は、点$(-2,0)$、$(-1,-a)$を通る

$W_1(s)$を
$\quad W_1(s)=As^3+Bs^2+Cs+D$
とすると、その1次微分は
$\quad W_1'(s)=3As^2+2Bs+C$
です。これらに値をあてはめると以下の4元連立方程式になります。
$\quad W_1(0)=D=1$
$\quad W_1(1)=A+B+C+D=0$
$\quad W_1'(0)=C=0$
$\quad W_1'(1)=3A+2B+C=a$
これを解くと
$\quad A=a+2$
$\quad B=-(a+3)$
$\quad C=0$
$\quad D=1$
となり、$W_1(s)$は
$\quad W_1(s)=(a+2)s^3-(a+3)s^2+1$
となります。
同様に$W_0(s)$は
$\quad W_0(s)=as^3-5as^2+8as-4a$
となります。

区間 $x_j \le x \le x_{j+1}$ の補間関数$g(x)$は、入力データ $f(x_{j-1})$、$f(x_j)$、$f(x_{j+1})$、$f(x_{j+2})$ に上記それぞれを乗じた和です。式で表すと；
$\quad g(x)=W_0(s+1)\ f(x_{j-1}) + W_1(s)\ f(x_j) + W_2(s-1)\ f(x_{j+1}) + W_3(s-2)\ f(x_{j+2})$

$W_2(s)$は$W_1(s)$の$s$の符号反転、$W_3(s)$は$W_0(s)$の$s$の符号反転なので；
$\quad g(x)=W_0(s+1)\ f(x_{j-1}) + W_1(s)\ f(x_j) + W_1(1-s)\ f(x_{j+1}) + W_0(2-s)\ f(x_{j+2})`$

任意の間隔の$x$について$g(x)$を求めることで再サンプリングとなります。画像処理ではこれを$x$方向と$y$方向に順次適用します。（順序は任意）
間隔が一定で$h$より小さいと画像拡大、大きいと画像縮小、等しいがオフセットしていると画像シフトになります。透視変換やモーフィングなどでは間隔は可変です。

行列表記

画像処理のバイキュービック補間については、畳み込みカーネル関数を変数変換すると；
$\quad W_0(s+1)=as^3-2as^2+as$
$\quad W_1(1-s)=-(a+2)s^3+(2a+3)s^2-as$
$\quad W_0(2-s)=-as^3+as^2$
これを用いて補間関数$g(x)$の式を行列表記にまとめると；

g(x)=
\begin{bmatrix}
s^3 & s^2 & s & 1
\end{bmatrix}
\begin{bmatrix}
a & a+2 & -a-2 & -a \\
-2a & -a-3 & 2a+3 & a \\
a & 0 & -a & 0 \\
0 & 1 & 0 & 0
\end{bmatrix}
\begin{bmatrix}
f(x_{j-1}) \\
f(x_j) \\
f(x_{j+1}) \\
f(x_{j+2})
\end{bmatrix}

Catmull-Romスプラインとの等価性

Catmull-Romスプラインはデザインやアニメーションのキーフレーム補間で用いられる曲線表現の一種と認知されていますが、数学的には画像処理のバイキュービック補間と等価です。

Catmull-Romスプラインは3次エルミートスプラインの一種です。3次エルミートスプラインでは、2点$\boldsymbol{p}_{k}$、$\boldsymbol{p}_{k+1}$の間が、接ベクトル$\boldsymbol{m}_{k}$、$\boldsymbol{m}_{k+1}$を用いる以下の式により補間されます。
$\quad \boldsymbol{p}(t) = \left(2t^3 - 3t^2 + 1\right) \boldsymbol{p}_{k} + \left(t^3 - 2t^2 + t\right) \boldsymbol{m}_{k} + \left(-2t^3 + 3t^2\right) \boldsymbol{p}_{k+1} + \left(t^3 - t^2\right) \boldsymbol{m}_{k+1}$
ただし、$t$はパラメトリック変数（$t \in [0,1]$）です。

接ベクトルをどう設定するかにより、さまざまな種類のスプラインを定義できます。接ベクトルを

\boldsymbol{m}_k = \frac{1}{2}(\boldsymbol{p}_{k+1} - \boldsymbol{p}_{k-1})

としたものがCatmull-Romスプラインです。
接ベクトル $\boldsymbol{m}_k = \frac{1}{2}(\boldsymbol{p}_{k+1} - \boldsymbol{p}_{k-1})$ を3次エルミートスプラインの式に代入して行列表記にまとめると；

\boldsymbol{p}(t) =
\frac{1}{2}
\begin{bmatrix}
t^3 & t^2 & t & 1
\end{bmatrix}
\begin{bmatrix}
-1 & 3 & -3 & 1 \\
2 & -5 & 4 & -1 \\
-1 & 0 & 1 & 0 \\
0 & 2 & 0 & 0
\end{bmatrix}
\begin{bmatrix}
\boldsymbol{p}_{k-1} \\
\boldsymbol{p}_k \\
\boldsymbol{p}_{k+1} \\
\boldsymbol{p}_{k+2}
\end{bmatrix}

この行列表記を画像処理のバイキュービック補間の行列表記と比較すると、$a=-0.5$ とした場合に一致することがわかります。

一部の文献では、$\frac{1}{2}$をテンションパラメータ$\tau$と係数化して；

\boldsymbol{m}_k = \tau (\boldsymbol{p}_{k+1} - \boldsymbol{p}_{k-1})

としている場合もあります。この場合の行列表記は；

\boldsymbol{p}(t) =
\begin{bmatrix}
t^3 & t^2 & t & 1
\end{bmatrix}
\begin{bmatrix}
-\tau & 2-\tau & \tau-2 & \tau \\
2\tau & \tau-3 & 3-2\tau & -\tau \\
-\tau & 0 & \tau & 0 \\
0 & 1 & 0 & 0
\end{bmatrix}
\begin{bmatrix}
\boldsymbol{p}_{k-1} \\
\boldsymbol{p}_k \\
\boldsymbol{p}_{k+1} \\
\boldsymbol{p}_{k+2}
\end{bmatrix}

この行列表記を画像処理のバイキュービック補間の行列表記と比較すると、$a=-\tau$ とした場合に一致することがわかります。

バイキュービック補間の1パス処理

画像処理のバイキュービック補間は、畳み込みカーネルをX方向とY方向に順次適用する2パス処理として実装することが多いと思います。一方で、Catmull-Romスプラインは、当初から二次元の補間（バイキュービックパッチ）に拡張できることが示されています。したがって、画像処理のバイキュービック補間をCatmull-Romスプラインのバイキュービックパッチとして実装すれば1パス処理が実現できますFPGAやASICなどのハード実装やGPUのような処理系では、このほうが効率の良い実装ができるのではないかと思います。

Catmull-Romスプラインのバイキュービックパッチの定義式は；

\boldsymbol{p}(s,t) =
\begin{bmatrix}
s^3 & s^2 & s & 1
\end{bmatrix}
\boldsymbol{M}
\begin{bmatrix}
\boldsymbol{p}_{11} & \boldsymbol{p}_{12} & \boldsymbol{p}_{13} & \boldsymbol{p}_{14} \\
\boldsymbol{p}_{21} & \boldsymbol{p}_{22} & \boldsymbol{p}_{23} & \boldsymbol{p}_{24} \\
\boldsymbol{p}_{31} & \boldsymbol{p}_{32} & \boldsymbol{p}_{33} & \boldsymbol{p}_{34} \\
\boldsymbol{p}_{41} & \boldsymbol{p}_{42} & \boldsymbol{p}_{43} & \boldsymbol{p}_{44}
\end{bmatrix}
\boldsymbol{M}^T
\begin{bmatrix}
t^3 \\
t^2 \\
t \\
1 \\
\end{bmatrix}

ただし

\boldsymbol{M}=
\frac{1}{2}
\begin{bmatrix}
-1 & 3 & -3 & 1 \\
2 & -5 & 4 & -1 \\
-1 & 0 & 1 & 0 \\
0 & 2 & 0 & 0
\end{bmatrix}

これを画像処理のバイキュービック補間に用いれば、1パス処理ができるはずです。

参考文献

Robert G. Keys
Cubic convolution interpolation for digital image processing
IEEE Transactions on Acoustics, Speech, and Signal Processing Vol. ASSP-29, No.6, 1981, pp.1153-1160
George Wolberg
Digital Image Warping
ISBN 0-8186-8944-7, pp.129-131
Edwin Catmull, Raphael Rom
A class of local interpolating splines
Computer Aided Geometric Design - Proceedings of a Conference Held at The University of Utah, Salt Lake City, Utah, March 18-21, 1974
ISBN 0-12-079050-5, pp.317-326
Wikipedia -- Cubic Hermite spline
https://en.wikipedia.org/wiki/Cubic_Hermite_spline
Christopher Twigg
Catmull-Rom splines
Carnegie Mellon Computer Graphics, 2003
https://www.cs.cmu.edu/~fp/courses/graphics/asst5/catmullRom.pdf

関連記事