【画像処理100本ノックに挑戦】Q.38. JPEG圧縮 (Step.2)DCT+量子化

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

# Q.38. JPEG圧縮 (Step.2)DCT+量子化
> DCT係数を量子化し、IDCTで復元せよ。また、その時の画像の容量を比べよ。

>Q = np.array(((16, 11, 10, 16, 24, 40, 51, 61),
              (12, 12, 14, 19, 26, 58, 60, 55),
              (14, 13, 16, 24, 40, 57, 69, 56),
              (14, 17, 22, 29, 51, 87, 80, 62),
              (18, 22, 37, 56, 68, 109, 103, 77),
              (24, 35, 55, 64, 81, 104, 113, 92),
              (49, 64, 78, 87, 103, 121, 120, 101),
              (72, 92, 95, 98, 112, 100, 103, 99)), dtype=np.float32)

解くことは難しくありません。係数が大きいと丸められやすいので、低周波ほど値が小さくなっています。上の数表をコードにするのが面倒ですが、vimの文字列置換で対応しました。
なお、画像の容量については計算していません。本来、符号化しないと容量は変わらないはずだからです。

```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 = 8;

	double q[64] = { 16, 11, 10, 16, 24, 40, 51, 61,
			  12, 12, 14, 19, 26, 58, 60, 55,
			  14, 13, 16, 24, 40, 57, 69, 56,
			  14, 17, 22, 29, 51, 87, 80, 62,
			  18, 22, 37, 56, 68, 109, 103, 77,
			  24, 35, 55, 64, 81, 104, 113, 92,
			  49, 64, 78, 87, 103, 121, 120, 101,
			  72, 92, 95, 98, 112, 100, 103, 99 };
	double Q[8][8];
	for (int j = 0; j < 8; j++)
		for (int i = 0; i < 8; i++)
		{
			Q[i][j] = q[i + 8 * j];
		}



	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);
				for (int i = 0; i < w; i++)
					for (int j = 0; j < h; j++)
					{
						out[i][j] = std::round(out[i][j] / Q[i][j]) * Q[i][j];
					}
				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/3a6cc31b-2dbf-d951-9f42-df9e75bbe15c.png) ![in.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520857/9db076b0-6fce-33e7-a5be-654ec5c6dc17.png)