はじめに

バイキュービック補間は、画像の画素と画素の間を、 3次関数によりXおよびY方向に補間するものです。画素と画素の間の値を求めることで、画像をずらす（シフトする）ことや、拡大することができます。

補間関数

各画素について補間計算をするのは大変なので、あらかじめ補間関数を求めておき、入力画像の各画素に畳み込むこと（コンボリューション）で補間した値を求めます。
バイキュービック補間の補間関数は3次のスプライン補間と似ていますが、2次導関数が連続との条件がないことが異なります。その分自由度があり、特性を調節することができます。

補間関数はデルタ関数画像（$x=0$だけ$1$、他は$0$）の画素と画素の間を補間する以下の3次関数の合成です。
- $x$が1～2の区間については、$(1,0)$,$(2,0)$を通る$f_0(x)$
- $x$が0～1の区間については、$(0,1)$,$(1,0)$を通る$f_1(x)$
- $x$が-1～0の区間については、$(-1,0)$,$(0,1)$を通る$f_2(x)$
- $x$が-2～-1の区間については、$(-2,0)$,$(-1,0)$を通る$f_3(x)$
補間関数の1次導関数が連続となるよう、以下の条件を追加します。（$A$は定数）
- $f_0(x)$の1次導関数$f_0'(x)$は$(1,A)$,$(2,0)$を通る
- $f_1(x)$の1次導関数$f_1'(x)$は$(0,0)$,$(1,A)$を通る
- $f_2(x)$の1次導関数$f_2'(x)$は$(-1,-A)$,$(0,0)$を通る
- $f_3(x)$の1次導関数$f_3'(x)$は$(-2,0)$,$(-1,-A)$を通る

$f_1(x)$を
　$f_1(x)=ax^3+bx^2+cx+d$
とすると、その1次導関数は
　$f_1'(x)=3ax^2+2bx+c$
です。これらに値をあてはめると以下の4元連立方程式になります。
　$f_1(0)=d=1$
　$f_1(1)=a+b+c+d=0$
　$f_1'(0)=c=0$
　$f_1'(1)=3a+2b+c=A$
これを解くと
　$a=A+2$
　$b=-(A+3)$
　$c=0$
　$d=1$
となり、$f_1(x)$は
　$f_1(x)=(A+2)x^3-(A+3)x^2+1$
となります。
同様に$f_0(x)$は
　$f_0(x)=Ax^3-5Ax^2+8Ax-4A$
となります。
$f_2(x)$は$f_1(x)$の$x$の符号反転、$f_3(x)$は$f_0(x)$の$x$の符号反転です。
グラフは以下の通りです。$A$の値により形が異なります。
　　
1次導関数は以下の通りです。
　$f_0'(x)=3Ax^2-10Ax+8A$
　$f_1'(x)=3(A+2)x^2-2(A+3)x$
　　
2次導関数は以下の通りです。$A=-0.75$の場合は連続になります。
　$f_0''(x)=6Ax-10A$
　$f_1''(x)=6(A+2)x-2(A+3)$

入力画像に補間関数を畳み込む

入力画像の近傍4画素の値に、補間関数を平行移動して乗じた合計が、補間された値です。

例として、$x$が$0$～$1$の区間については、$x=-1$の画素の値に$f_0(x+1)$を、$x=0$の画素の値に$f_1(x)$を、$x=1$の画素の値に$f_2(x-1)$を、$x=2$の画素の値に$f_3(x-2)$を、それぞれ乗じた合計が補間された値になります。

以下は任意の$x$について、補間された値を求める処理の例です。
（数式だとわかりにくくなるので、プログラムです）

double cubicInterpolatedValue(
  double *img,	// 入力画像（1次元）
  int xi,		// xの整数部
  double xf)	// xの小数部
{
	double rv;
	rv  = img[xi – 1] * f0(xf + 1.0);
	rv += img[xi]     * f1(xf);
	rv += img[xi + 1] * f2(xf – 1.0);
	rv += img[xi + 2] * f3(xf – 2.0);
	return rv;
}

2次元に拡張

X方向に補間し、次にY方向に補間することで2次元の画像補間ができます。
補間関数を2次元に拡張して、X方向とY方向をいちどに補間することもできます。
2次元に拡張した補間関数は以下のようになります。

実装例：OpenCVの補間関数

OpenCVの補間関数は以下の通りです。（4.6.0からの抜粋。他バージョンも同一です）

modules/imgproc/src/resize.cpp から抜粋

static inline void interpolateCubic( float x, float* coeffs )
{
    const float A = -0.75f;

    coeffs[0] = ((A*(x + 1) - 5*A)*(x + 1) + 8*A)*(x + 1) - 4*A;
    coeffs[1] = ((A + 2)*x - (A + 3))*x*x + 1;
    coeffs[2] = ((A + 2)*(1 - x) - (A + 3))*(1 - x)*(1 - x) + 1;
    coeffs[3] = 1.f - coeffs[0] - coeffs[1] - coeffs[2];
}

加算や減算が括弧でくくられているのは計算誤差を抑えるためです。
（次数が大きく異なる項どうし足したり引いたりしない）
$f_3(x)$は$f_0(x)$～$f_3(x)$の和が1であること利用して計算量が低減されています。
$A$は$-0.75$固定です。

※OpenCVのビルドの設定によっては、他のコードが用いられる場合があります。
※OpenCVのライセンス条件は以下の通りです。

参考文献

George Wolberg
Digital Image Warping, pp.129-131
ISBN 0-8186-8944-7

関連記事