バイキュービック補間

数学において、**バイキュービック補間**（バイキュービックほかん、双3次補間、）は2次元の上のデータポイントを補間するための3次スプライン補間（データセットに3次補間を適用する方法）の拡張である。
補間された曲面 (画像ではなくカーネルの形状) は、バイリニア補間または最近傍補間によって得られる曲面よりも滑らかになる。
バイキュービック補間は、ラグランジュ多項式、3次スプライン補間、または3次畳み込みアルゴリズムを使用して実現できる。

画像処理では、処理速度が問題にならない場合の画像のリサンプリングで、バイリニア補間や最近傍補間よりもバイキュービック補間が選択されることが多い。
4ピクセル（2×2）のみを考慮するバイリニア補間とは対照的に、バイキュービック補間では16ピクセル（4×4）を考慮する。
バイキュービック補間で再サンプリングされた画像は、選択されたb値とc値に応じて、さまざまな補間アーティファクトが発生する可能性がある。

## 計算

四隅が (0,0)、(1,0)、(0,1)、(1,1) の単位正方形の関数を f、
微分を f_x、f_y、f_{xy} とする。
補間された曲面は以下で表される。
p(x,y) = \sum\limits_{i=0}^3 \sum_{j=0}^3 a_{ij} x^i y^j

補間問題は16個の係数 a_{ij} を決定することから成り立つ。
p(x,y) に関数をあてはめ、4つの方程式を得る。
# f(0,0)      = p(0,0)   = a_{00}
# f(1,0)      = p(1,0)   = a_{00} + a_{10} + a_{20} + a_{30}
# f(0,1)      = p(0,1)   = a_{00} + a_{01} + a_{02} + a_{03}
# f(1,1)      = p(1,1)   = \textstyle \sum\limits_{i=0}^3 \sum\limits_{j=0}^3 a_{ij}

同様に、xおよびy方向の微分から、8つの方程式を得る。
# f_x(0,0)    = p_x(0,0) = a_{10}
# f_x(1,0)    = p_x(1,0) =  a_{10} + 2a_{20} + 3a_{30}
# f_x(0,1)    = p_x(0,1) = a_{10} + a_{11} + a_{12} + a_{13}
# f_x(1,1)    = p_x(1,1) = \textstyle \sum\limits_{i=1}^3 \sum\limits_{j=0}^3 a_{ij} i
# f_y(0,0)    = p_y(0,0) = a_{01}
# f_y(1,0)    = p_y(1,0) = a_{01} + a_{11} + a_{21} + a_{31}
# f_y(0,1)    = p_y(0,1) = a_{01} + 2a_{02} + 3a_{03}
# f_y(1,1)    = p_y(1,1) = \textstyle \sum\limits_{i=0}^3 \sum\limits_{j=1}^3 a_{ij} j

さらに、xyの混合偏微分から4つの方程式を得る。
# f_{xy}(0,0) = p_{xy}(0,0) = a_{11}
# f_{xy}(1,0) = p_{xy}(1,0) = a_{11} + 2a_{21} + 3a_{31}
# f_{xy}(0,1) = p_{xy}(0,1) = a_{11} + 2a_{12} + 3a_{13}
# f_{xy}(1,1) = p_{xy}(1,1) = \textstyle \sum\limits_{i=1}^3 \sum\limits_{j=1}^3 a_{ij} i j

上記表現で、次の恒等式が使用されている。
p_x(x,y) = \textstyle \sum\limits_{i=1}^3 \sum\limits_{j=0}^3 a_{ij} i x^{i-1} y^j
p_y(x,y) = \textstyle \sum\limits_{i=0}^3 \sum\limits_{j=1}^3 a_{ij} x^i j y^{j-1}
p_{xy}(x,y) = \textstyle \sum\limits_{i=1}^3 \sum\limits_{j=1}^3 a_{ij} i x^{i-1} j y^{j-1}

この手順により求められる単位正方形 [0,1] \times [0,1] 上の曲面 p(x,y) は連続、かつ導関数連続になる。
任意のサイズの規則格子上のバイキュービック補間は、このようなバイキュービック曲面をパッチし、境界上で導関数が一致するようにすることで実現できる。

未知のパラメータ a_{ij} をグループ化したベクトルを
\alpha=\left[\begin{smallmatrix}a_{00}&a_{10}&a_{20}&a_{30}&a_{01}&a_{11}&a_{21}&a_{31}&a_{02}&a_{12}&a_{22}&a_{32}&a_{03}&a_{13}&a_{23}&a_{33}\end{smallmatrix}\right]^T
とし、
x=\left[\begin{smallmatrix}f(0,0)&f(1,0)&f(0,1)&f(1,1)&f_x(0,0)&f_x(1,0)&f_x(0,1)&f_x(1,1)&f_y(0,0)&f_y(1,0)&f_y(0,1)&f_y(1,1)&f_{xy}(0,0)&f_{xy}(1,0)&f_{xy}(0,1)&f_{xy}(1,1)\end{smallmatrix}\right]^T
とすると、上記連立方程式は、線形方程式 A\alpha=x の行列に再定式化できる。

逆行列により、有用な線形方程式  A^{-1}x=\alpha が得られる。ここで、
A^{-1}=\left[\begin{smallmatrix}\begin{array}{rrrrrrrrrrrrrrrr}
 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 -3 & 3 & 0 & 0 & -2 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 2 & -2 & 0 & 0 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -3 & 3 & 0 & 0 & -2 & -1 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2 & -2 & 0 & 0 & 1 & 1 & 0 & 0 \\
 -3 & 0 & 3 & 0 & 0 & 0 & 0 & 0 & -2 & 0 & -1 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & -3 & 0 & 3 & 0 & 0 & 0 & 0 & 0 & -2 & 0 & -1 & 0 \\
 9 & -9 & -9 & 9 & 6 & 3 & -6 & -3 & 6 & -6 & 3 & -3 & 4 & 2 & 2 & 1 \\
 -6 & 6 & 6 & -6 & -3 & -3 & 3 & 3 & -4 & 4 & -2 & 2 & -2 & -2 & -1 & -1 \\
 2 & 0 & -2 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 2 & 0 & -2 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 \\
 -6 & 6 & 6 & -6 & -4 & -2 & 4 & 2 & -3 & 3 & -3 & 3 & -2 & -1 & -2 & -1 \\
 4 & -4 & -4 & 4 & 2 & 2 & -2 & -2 & 2 & -2 & 2 & -2 & 1 & 1 & 1 & 1
\end{array}\end{smallmatrix}\right]
である。これにより \alpha を素早く簡単に計算できる。

16個の係数に対し、別の簡潔な行列形式が存在する。
\begin{bmatrix}
f(0,0)&f(0,1)&f_y (0,0)&f_y (0,1)\\f(1,0)&f(1,1)&f_y (1,0)&f_y (1,1)\\f_x (0,0)&f_x (0,1)&f_{xy} (0,0)&f_{xy} (0,1)\\f_x (1,0)&f_x (1,1)&f_{xy} (1,0)&f_{xy} (1,1)
\end{bmatrix} =
\begin{bmatrix}
1&0&0&0\\1&1&1&1\\0&1&0&0\\0&1&2&3
\end{bmatrix}
\begin{bmatrix}
a_{00}&a_{01}&a_{02}&a_{03}\\a_{10}&a_{11}&a_{12}&a_{13}\\a_{20}&a_{21}&a_{22}&a_{23}\\a_{30}&a_{31}&a_{32}&a_{33}
\end{bmatrix}
\begin{bmatrix}
1&1&0&0\\0&1&1&1\\0&1&0&2\\0&1&0&3
\end{bmatrix}
または

\begin{bmatrix}
a_{00}&a_{01}&a_{02}&a_{03}\\a_{10}&a_{11}&a_{12}&a_{13}\\a_{20}&a_{21}&a_{22}&a_{23}\\a_{30}&a_{31}&a_{32}&a_{33}
\end{bmatrix}
=
\begin{bmatrix}
1&0&0&0\\0&0&1&0\\-3&3&-2&-1\\2&-2&1&1
\end{bmatrix}
\begin{bmatrix}
f(0,0)&f(0,1)&f_y (0,0)&f_y (0,1)\\f(1,0)&f(1,1)&f_y (1,0)&f_y (1,1)\\f_x (0,0)&f_x (0,1)&f_{xy} (0,0)&f_{xy} (0,1)\\f_x (1,0)&f_x (1,1)&f_{xy} (1,0)&f_{xy} (1,1)
\end{bmatrix}
\begin{bmatrix}
1&0&-3&2\\0&0&3&-2\\0&1&-2&1\\0&0&-1&1
\end{bmatrix},

ただし
p(x,y)=\begin{bmatrix}1
&x&x^2&x^3\end{bmatrix}
\begin{bmatrix}
a_{00}&a_{01}&a_{02}&a_{03}\\a_{10}&a_{11}&a_{12}&a_{13}\\a_{20}&a_{21}&a_{22}&a_{23}\\a_{30}&a_{31}&a_{32}&a_{33}
\end{bmatrix}
\begin{bmatrix}1\\y\\y^2\\y^3\end{bmatrix}
である。

## 不規則格子への拡張
多くのアプリケーションでは、単位正方形ではなく上のデータを使用してバイキュービック補間を行う。
この場合 p_x、p_y、p_{xy} の恒等式は

p_x(x,y) = \textstyle \sum\limits_{i=1}^3 \sum\limits_{j=0}^3 \frac{a_{ij} i x^{i-1} y^j}{\Delta x}
p_y(x,y) = \textstyle \sum\limits_{i=0}^3 \sum\limits_{j=1}^3 \frac{a_{ij} x^i j y^{j-1}}{\Delta y}
p_{xy}(x,y) = \textstyle \sum\limits_{i=1}^3 \sum\limits_{j=1}^3 \frac{a_{ij} i x^{i-1} j y^{j-1}}{\Delta x \Delta y}
となる。
ここで \Delta x と \Delta y は 点 (x,y) を含むセルの間隔である。
この場合、係数 \alpha を計算する最も実用的な方法は、

x=\left[\begin{smallmatrix}f(0,0)&f(1,0)&f(0,1)&f(1,1)&\Delta x f_x(0,0)&\Delta xf_x(1,0)&\Delta x f_x(0,1)&\Delta x f_x(1,1)&\Delta y f_y(0,0)&\Delta y f_y(1,0)&\Delta y f_y(0,1)&\Delta y f_y(1,1)&\Delta x \Delta y f_{xy}(0,0)&\Delta x \Delta y f_{xy}(1,0)&\Delta x \Delta y f_{xy}(0,1)&\Delta x \Delta y f_{xy}(1,1)\end{smallmatrix}\right]^T
として、前述の A を用いて \alpha=A^{-1}x を解くことである。次に正規化された補間変数を次のように計算する。

\begin{align}
\overline{x} &= \frac{x-x_0}{x_1-x_0}, \\
\overline{y} &= \frac{y-y_0}{y_1-y_0}
\end{align}
ここで x_0, x_1, y_0 と y_1 は点 (x,y) を囲むグリッド点の x および y 座標である。すると、補間された曲面は
p(x,y) = \sum\limits_{i=0}^3 \sum_{j=0}^3 a_{ij} {\overline{x}}^i {\overline{y}}^j
となる。

## 関数値から導関数を求める

導関数は、不明な場合は通常、有限差分を用いるなどして単位正方形の角に隣接する点の関数から近似される。

1次導関数 f_x または f_y のいずれかを求めるには、適切な軸上の*囲っている*2点間の傾きを計算する方法を用いる。
たとえば、ある1点の f_x を計算するには、目標点の左右の点について f(x,y) を求めて傾きを計算する。f_y についても同様にする。

交差導関数 f_{xy} を求めるには、両方の軸で一度に1つずつ導関数を得る。
例えば、まず f_x の手順を用いて、目標点の上と下の点の x の導関数を算出する。
次に、それら値に f_y の手順を用いて（通常のように、それらの点の f の値ではなく）目標点の f_{xy}(x,y) の値を計算する 。
（あるいは逆の順序で、最初に f_y を計算し、次にそれらから f_x を計算することもできる。両者とも同じ結果になる。）

データセットの端部で周囲のポイントの一部が欠落している場合、欠落しているポイントはいくつかの方法で近似できる。
単純で一般的な方法は、既存のポイントからターゲットポイントまでの傾斜が、それ以降変化せずに続くと仮定し、これを使用して欠落しているポイントの仮想値を計算することである。

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

バイキュービックスプライン補間では、各グリッドセルに対して上記の線形システムの解が必要である。
同様の特性を持つ補間は、両次元に次のカーネルを使用した畳み込みを適用することで求めることができる。

W(x) = 
\begin{cases}
 (a+2)|x|^3-(a+3)|x|^2+1 & \text{for } |x| \leq 1 \\
 a|x|^3-5a|x|^2+8a|x|-4a & \text{for } 1 
ここで a は通常、−0.5または−0.75に設定される。
なお、W(0)=1、かつ、すべての0以外の整数 n に対し W(n)=0 である。

このアプローチは ERTS（地球資源技術衛星、のちのランドサット1号）の画像データを対象に、TRW Systems GroupのRifmanにより提案された

。
当初 a は-1固定であったが、同社Simonによりパラメータ化され、-0.5、-0.75、-1.0が意味のある値とされた
。
しかしこれら提案では導出過程や数式などが十分に提示されていなかったため、のちに Keys によって完全な形で再提案され、a=-0.5 の場合に、元の関数のサンプリング間隔について3次収束が得られることが示された
。

畳み込みカーネル W(x) は以下により導出される。
W(x) = 
\begin{cases}
 A_1|x|^3+B_1|x|^2+C_1|x|+D_1 & \text{for } |x| \leq 1 \\
 A_2|x|^3+B_2|x|^2+C_2|x|+D_2 & \text{for } 1 
とする。W(0)=1、W(1)=W(2)=0、また1次微分については W'(0)=W'(2)=0、W'(1)=a なので、
# W(0)=D_1=1
# W(1^-)=A_1+B_1+C_1+D_1=0
# W(1^+)=A_2+B_2+C_2+D_2=0
# W(2^-)=8A_2+4B_2+2C_2+D_2=0
# W'(0)=C_1=0
# W'(1^-)=3A_1+2B_1+C_1=a
# W'(1^+)=3A_2+2B_2+C_2=a
# W'(2^-)=12A_2+4B_2+C_2=0

以上の連立方程式を解くと、
A_1=a+2, \quad B_1=-(a+3), \quad C_1=0, \quad D_1=1
A_2=a, \quad B_2=-5a, \quad C_2=8a, \quad D_2=-4a
となり、上記式を得る。

### 行列表記

行列表記を使用すると、方程式をよりわかりやすい方法で表現できる。
ひとつの次元で0から1の間の t について、
p(t) =
\begin{bmatrix}
1 & t & t^2 & t^3
\end{bmatrix}

\begin{bmatrix}
0 & 1 & 0 & 0 \\
a & 0 & -a & 0 \\
-2a & -a-3 & 2a+3 & a \\
a & a+2 & -a-2 & -a
\end{bmatrix}

\begin{bmatrix}
f_{-1} \\
f_0 \\
f_1 \\
f_2
\end{bmatrix}

となる。
1次元の3次畳み込み補間では4つのサンプル点が必要である。
計算ごとに、その左側の2つのサンプルと、右側の2つのサンプルが参照される。
このテキストでは、これらのポイントに-1から2までインデックス付けされている。
ここでは、0でインデックス付けされたポイントから照会ポイントまでの距離は t で示される。

a=-0.5 の場合、Catmull-Romスプラインと一致する。
このため、a を -0.5 に設定したバイキュービック補間を "Catmull-Rom interpolation" と呼ぶ場合がある
。
また、a=-0.5 以外の場合は、テンションパラメータが拡張されたCatmull-Romスプラインと一致し、テンションパラメータ \tau=-a である。

2つの次元について、最初に x、次に y に適用する場合、
\begin{align}
b_{-1}&= p(t_x, f_{(-1,-1)}, f_{(0,-1)}, f_{(1,-1)}, f_{(2,-1)}), \\[1ex]
b_{0} &= p(t_x, f_{(-1,0)}, f_{(0,0)}, f_{(1,0)}, f_{(2,0)}), \\[1ex]
b_{1} &= p(t_x, f_{(-1,1)}, f_{(0,1)}, f_{(1,1)}, f_{(2,1)}), \\[1ex]
b_{2} &= p(t_x, f_{(-1,2)}, f_{(0,2)}, f_{(1,2)}, f_{(2,2)}),
\end{align}
p(x,y) = p(t_y, b_{-1}, b_{0}, b_{1}, b_{2})
である。

### 微分連続性
定義により1次微分連続であるが、行列表記を用いた以下の方法でも容易に判定できる。
:
p_i'(1) =
\begin{bmatrix}
0 & a & 0 & -a
\end{bmatrix}
\begin{bmatrix}
f_{i-1} \\
f_i \\
f_{i+1} \\
f_{i+2}
\end{bmatrix}
=
p_{i+1}'(0) =
\begin{bmatrix}
a & 0 & -a & 0
\end{bmatrix}
\begin{bmatrix}
f_i \\
f_{i+1} \\
f_{i+2} \\
f_{i+3}
\end{bmatrix}

したがい、1次微分連続である。
:
p_i*(1) =
\begin{bmatrix}
2a & 4a+6 & -2a-6 & -4a
\end{bmatrix}
\begin{bmatrix}
f_{i-1} \\
f_i \\
f_{i+1} \\
f_{i+2}
\end{bmatrix}
\neq
p_{i+1}*(0) =
\begin{bmatrix}
-4a & -2a-6 & 4a+6 & 2a
\end{bmatrix}
\begin{bmatrix}
f_i \\
f_{i+1} \\
f_{i+2} \\
f_{i+3}
\end{bmatrix}

したがい、2次微分連続ではない。

### 他の方式との比較
以下は各方式のカーネル関数の比較である。

## コンピュータグラフィックスでの使用

バイキュービック アルゴリズムは、画像やビデオを表示用に拡大するためによく使用される (画像スケーリングを参照)。
一般的なバイリニアアルゴリズムよりも細かいディテールが保持される。

バイキュービック アルゴリズムは、画像縮小にも使用されるが、縮小時の再サンプリングは実光学系のズームアウトの作用と異なるため、他の方式のほうが適している場合がある。

### パラメータaの影響

畳み込みカーネルの負のローブにより、オーバーシュート (電気工学)（ハロー）が発生する。
これによりクリッピングが発生する可能性があり、アーティファクト（リンギングアーティファクトも参照）になるが、尖鋭度（見かけの鮮明度）が増すため、望ましい場合もある。この効果は、パラメータaが-0.5の場合より、-0.75や-1.0の場合のほうが大きい。

a=-0.5 の場合、補間結果が、元の画像（サンプリング前の連続関数）の2次微分までのテイラー近似（剰余項除く）と一致する
。
この作用を端的に示す例が、元の画像が1次関数（単純なグラデーションパターン）の場合であり、a=-0.5 の場合のみ、補間結果が1次関数になる。
例えば、入力が
\begin{bmatrix} f_{i-1} & f_i & f_{i+1} & f_{i+2} \end{bmatrix} = \begin{bmatrix} -1 & 0 & 1 & 2 \end{bmatrix}
であると、補間結果は
p(t) = -2(2a+1)t^3+3(2a+1)t^2-2at
であり、a=-0.5 の場合のみ1次関数（2次と3次の項が0）になる。

以下は単純なグラデーション画像の16倍拡大と、標準テスト画像"cameraman"の8倍拡大（一部）である。
前者の縦縞状の濃淡むらや、後者のこめかみ部分の人工的なまだら模様は、上記作用によるアーティファクトである。

なお、a=-0.75 の場合は、畳み込みカーネルの2次微分が x=1 で連続になるが、補間結果が2次微分連続になるわけではない。
また、a=-1.0 の場合は、畳み込みカーネルの微分が x=1 でsinc関数の微分と一致するが、sinc補間との比較以外の意義はない。

### 再サンプリング位置
整数倍の画像拡大では、拡大画像の重心が入力画像に対しずれないよう、入力画像のピクセルをはさむ位置に再サンプリング位置を設定する（下記図左）ことが一般的である。
一方で、入力画像のピクセルに再サンプリング位置を重ねるよう設定する方法もある（下記図右）。この方法では変換の可逆性が確保される。

### 外挿補間
補間のために左右（または上下）に入力画像が2ピクセルずつ必要であり、画像の外周部は本来は無効な領域であるため、外挿補間の対応が必要になる。外挿補間の方式は実装により異なる。
以下は、最外周のピクセルをコピーして外挿する場合の例である。

## 脚注

## 関連項目

* 
* ベジェ曲面
* 画像スケーリング
* 最近傍補間
* バイリニア補間
* 
* 
* 
* 
* 
* 
* 

## 外部リンク
* 
* Interpolation theory
* Explanation and Java/C++ implementation of (bi)cubic interpolation
* Excel Worksheet Function for Bicubic Lagrange Interpolation