使用したライブラリ
Q.27. Bi-cubic補間
Bi-cubic補間により画像を1.5倍に拡大せよ。
前回の記事(サンプリング定理と画像処理における補間)にて、離散信号から元の連続信号を復元するときにはsinc関数を畳み込めばよいことを確認しました。バイキュービック補間はsinc関数の近似式を用いているらしいですが本当でしょうか。確かめてみます。
最初はテイラー展開してやればよいかとも思ったのですが、sinc関数のテイラー展開は偶数次しか持ちません。どうもそこまで話は単純でないようです。調べてみるとどうやら、$0\leq x\leq 1,1< x< 2$の範囲でそれぞれ異なる3次関数を用いるらしいです。つまり
\phi(x) =
\begin{cases}
\sum_{i=0}^3 a_ix^i & 0\leq x\leq 1\\
\sum_{i=0}^3 b_ix^i & 1< x< 2\\
0 & {\rm otherwise}
\end{cases}
とします。対称な形を考えているので、$x$が負の領域も同様に定義されます。sinc関数の近似とするために、$\phi(0)=1,\phi(1)=\phi(2)=0$が必要です。さらに、領域を3つに分けているので、それぞれが滑らかに接続する必要があります。これらの条件から7つの連立方程式を立てられますが、ひとつだけ任意パラメータが残ります。それが問題文の式になり、パラメータは$a$と置かれています。多くの場合$a=-1$が用いられるそうですが、本当にsinc関数の近似になっているか確認します。
何をもって良い近似とするかは主観が入りますが、ここでは簡単に二乗誤差を指標としてみます。
\int_0^2 \{\phi(x;a)-\mathrm{sinc}x\}^2dx
解析的に解けるのかもしれませんが、面倒なので数値的にごり押しします。
# include <iostream>
# include <fstream>
double pi = 3.14159265358979;
int main()
{
auto sinc = [](double x)
{
if (x < 1e-5) return 1.;
return sin(pi * x) / pi / x;
};
auto phi = [](double x,double a)
{
if (0 <= x && x < 1) return (a + 2) * x * x * x - (a + 3) * x * x + 1;
else if (1 <= x && x < 2) return a * x * x * x - 5 * a * x * x + 8 * a * x - 4 * a;
else return 0.;
};
double amin;
double evalmin = 1e15;
for (double a = -5; a <= 5; a += 0.001)
{
double sum = 0;
for (double x = 0; x <= 2; x += 0.0001)
{
double val = (sinc(x) - phi(x, a));
sum += val * val;
}
if (sum < evalmin)
{
evalmin = sum;
amin = a;
}
std::cout << a << " " << sum << std::endl;
}
std::cout << amin << std::endl;
std::ofstream ofs0("-1.txt");
std::ofstream ofs("amin.txt");
for (double x = 0; x < 3; x += 0.01)
{
ofs0 << x << " " << phi(x, -1) << std::endl;
ofs << x << " " << phi(x, amin) << std::endl;
}
}
二乗誤差的には$a=-1.295$の方がよさそうです。でもそれほど大きくは変わりませんね。$a=-1$もそれなりにsinc関数の近似として使えそうです。でもいっそsinc関数そのものを用いた方が早くない？と思うのですが。しかし、sinc関数をx=2で切ってしまうと傾きが滑らかに接続しませんね。そこが問題かもしれません。 Lanczos補間というものもあるようですが、これはsinc関数に窓関数をかけているので、傾きは滑らかなようです。
それではバイキュービック補間をやってみます。
int main()
{
PPM ppm("imori.pnm");
int width = ppm.Get_width();
int height = ppm.Get_height();
double a = 1.5;
int Width = width * a, Height = height * a;
auto phi = [](double t)
{
double x = fabs(t);
double a = -1;
if (0 <= x && x < 1) return (a + 2) * x * x * x - (a + 3) * x * x + 1;
else if (1 <= x && x < 2) return a * x * x * x - 5 * a * x * x + 8 * a * x - 4 * a;
else return 0.;
};
PPM ppm2(Width, Height);
for(int J=0; J<Height; J++)
for (int I = 0; I < Width; I++)
{
ppm2(I, J, 'r') = 0;
ppm2(I, J, 'g') = 0;
ppm2(I, J, 'b') = 0;
double x = I / a, y = J / a;
int i = x, j = y;
double norm = 0;
double r, g, b;
r = g = b = 0;
for(int dj=-2; dj<=2; dj++)
for (int di = -2; di <= 2; di++)
{
double tx = fabs(i + di - x);
double ty = fabs(j + dj - y);
if ((i + di < 0 || i + di >= width) || (j + dj < 0 || j + dj >= height))
{
}
else
{
norm += phi(tx) * phi(ty);
r += ppm(i + di, j + dj, 'r') * phi(tx) * phi(ty);
g += ppm(i + di, j + dj, 'g') * phi(tx) * phi(ty);
b += ppm(i + di, j + dj, 'b') * phi(tx) * phi(ty);
}
}
ppm2(I, J, 'r') = (r / norm <=255 ? r/norm : 255);
ppm2(I, J, 'g') = (g / norm <= 255 ? g / norm : 255);
ppm2(I, J, 'b') = (b / norm <= 255 ? b / norm : 255);
if (r/norm > 255) std::cout << norm << std::endl;
}
ppm2.Flush("out.ppm");
return 0;
}
境界の扱い等々いい加減なので注意。ちなみに$a=-1.295$でもやってみましたがほぼ同じでした。各色8bitしかないですしね。