【画像処理100本ノックに挑戦】Q.27. Bi-cubic補間

# 使用したライブラリ
[【画像処理100本ノック】独自の画像入出力クラスを作る](https://qiita.com/jajagacchi/items/834e8425b84f94df8a49)

# Q.27. Bi-cubic補間
> Bi-cubic補間により画像を1.5倍に拡大せよ。

前回の記事([サンプリング定理と画像処理における補間](https://qiita.com/jajagacchi/items/1d1fcc852ef7d470fe7e))にて、離散信号から元の連続信号を復元するときにはsinc関数を畳み込めばよいことを確認しました。バイキュービック補間はsinc関数の近似式を用いているらしいですが本当でしょうか。確かめてみます。

最初はテイラー展開してやればよいかとも思ったのですが、sinc関数のテイラー展開は偶数次しか持ちません。どうもそこまで話は単純でないようです。調べてみるとどうやら、$0\leq x\leq 1,1< x< 2$の範囲でそれぞれ異なる3次関数を用いるらしいです。つまり

```math
\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関数の近似になっているか確認します。

何をもって良い近似とするかは主観が入りますが、ここでは簡単に二乗誤差を指標としてみます。

```math
\int_0^2 \{\phi(x;a)-\mathrm{sinc}x\}^2dx
```

解析的に解けるのかもしれませんが、面倒なので数値的にごり押しします。

```c++
#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;
	}


}
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520857/b9a7dc0b-83d5-666f-a436-cbbbe8bcee0b.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520857/aeedc443-865e-fb01-0ac7-7dac4ac6f7ee.png)

二乗誤差的には$a=-1.295$の方がよさそうです。でもそれほど大きくは変わりませんね。$a=-1$もそれなりにsinc関数の近似として使えそうです。でもいっそsinc関数そのものを用いた方が早くない？と思うのですが。しかし、sinc関数をx=2で切ってしまうと傾きが滑らかに接続しませんね。そこが問題かもしれません。 Lanczos補間というものもあるようですが、これはsinc関数に窓関数をかけているので、傾きは滑らかなようです。

それではバイキュービック補間をやってみます。

```c++
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;
}
```
![imori.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520857/1fc980a9-c73b-1f3e-294f-9df45ef963c2.png) ![out.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520857/c8cc1d1b-e775-4967-ac72-c22211c04c7d.png)

境界の扱い等々いい加減なので注意。ちなみに$a=-1.295$でもやってみましたがほぼ同じでした。各色8bitしかないですしね。