【画像処理100本ノックに挑戦】Q.37. PSNR

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

# Q.37. PSNR
> IDCTで用いるDCT係数を8でなく、4にすると画像の劣化が生じる。 入力画像とIDCT画像のPSNRを求めよ。また、IDCTによるビットレートを求めよ。

簡単ですね。

```c++
std::vector<std::vector<double>> dct2d(const std::vector<std::vector<double>>& x, int M, int N)
{
	std::vector<std::vector<double>> ret(M, std::vector<double>(N));
	for (int m = 0; m < M; m++)
		for (int n = 0; n < N; n++)
		{
			ret[n][m] = 0;
		}
	for (int k = 0; k < M; k++)
		for (int l = 0; l < N; l++)
		{
			for (int m = 0; m < M; m++)
				for (int n = 0; n < N; n++)
				{
					ret[k][l] += x[m][n] * cos(pi / M * (m + 0.5) * k) * cos(pi / N * (n + 0.5) * l);
				}
		}
	return ret;
}

std::vector<std::vector<double>> idct2d(const std::vector<std::vector<double>>& x, int M, int N, int K)
{
	std::vector<std::vector<double>> ret(M, std::vector<double>(N));
	for (int m = 0; m < M; m++)
		for (int n = 0; n < N; n++)
		{
			ret[n][m] = 0;
		}
	for (int k = 0; k < M; k++)
		for (int l = 0; l < N; l++)
		{
			double a, b;
			for (int m = 0; m < K; m++)
				for (int n = 0; n < K; n++)
				{
					if (m == 0) a = 0.5;
					else a = cos(pi / M * m * (k + 0.5));
					if (n == 0) b = 0.5;
					else b = cos(pi / N * n * (l + 0.5));

					ret[k][l] += x[m][n] * a * b / (M * N / 4.);
				}
		}
	return ret;
}

int main()
{
	PPM ppm("imori.pnm");
	int M = ppm.Get_width();
	int N = ppm.Get_height();

	char c[3] = { 'r','g','b' };
	PPM ppm2(M, N);
	int K = 4;
	for (int color = 0; color < 3; color++)
	{
		std::vector<std::vector<double>> f(M, std::vector<double>(N));
		for (int j = 0; j < N; j++)
			for (int i = 0; i < M; i++)
			{
				f[i][j] = ppm(i, j, c[color]);
			}

		std::vector<std::vector<double>> f2(M, std::vector<double>(N));
		int w = 8, h = 8;
		
		for (int I = 0; I < M / w; I++)
			for (int J = 0; J < N / h; J++)
			{
				std::vector<std::vector<double>> in(w, std::vector<double>(h));
				for (int i = w * I; i < w * I + w; i++)
					for (int j = h * J; j < h * J + h; j++)
					{
						in[i - w * I][j - h * J] = f[i][j];
					}
				auto out = dct2d(in, w, h);
				auto in2 = idct2d(out, w, h, K);
				for (int i = w * I; i < w * I + w; i++)
					for (int j = h * J; j < h * J + h; j++)
					{
						ppm2(i, j, c[color]) = in2[i - w * I][j - h * J];
						if (ppm2(i, j, c[color]) < 0) ppm2(i, j, c[color]) = 0;
						else if (ppm2(i, j, c[color]) > 255) ppm2(i, j, c[color]) = 255;
					}
			}
	}

	double mse = 0;
	int maxval = 0;
	for(int i=0; i<M; i++)
		for (int j = 0; j < N; j++)
		{
			if (maxval < ppm2(i, j, 'r')) maxval = ppm2(i, j, 'r');
			if (maxval < ppm2(i, j, 'g')) maxval = ppm2(i, j, 'g');
			if (maxval < ppm2(i, j, 'b')) maxval = ppm2(i, j, 'b');
			mse += (ppm2(i, j, 'r') - ppm(i, j, 'r')) * (ppm2(i, j, 'r') - ppm(i, j, 'r')) / (double)(M * N * 3);
			mse += (ppm2(i, j, 'g') - ppm(i, j, 'g')) * (ppm2(i, j, 'g') - ppm(i, j, 'g')) / (double)(M * N * 3);
			mse += (ppm2(i, j, 'b') - ppm(i, j, 'b')) * (ppm2(i, j, 'b') - ppm(i, j, 'b')) / (double)(M * N * 3);
		}
	double psnr = 10 * log10(maxval * maxval / mse);
	double bitrate = 8 * K * K / 8./8.;
	std::cout << psnr << " " << bitrate;

	ppm2.Flush("in.pnm");
	return 0;
}
```
![imori.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520857/fcda9a77-cb54-bb9a-ee94-95d568e3ab87.png) ![in.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520857/121bf3f1-4699-4d3c-4a4e-310169173861.png)

8x8画素ごとにでDCTしているので、どうもこれに由来するアーチファクトがありますね(8x8のパッチが見える)。正解画像にもあるのでこんなもんみたいです。
少し調べてみましたが8x8画素ごとに処理を行うはっきりとした理由はよくわかりませんでした。計算時間はひとつの理由かもしれません。128x128でもDCTしてみたのですが、えらく計算時間が違います。上のコードの場合、大雑把に

```math
w^4\times\left(\frac{M}{w}\right)^2
```

の計算回数が必要です。$M=128$として$w=128$の場合と$w=8$の場合を比べると、後者は前者に比べて1/256の計算量で済みます。