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

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

# Q.40. JPEG圧縮 (Step.4)YCbCr+DCT+量子化
> YCbCr表色系にし、DCT後、Yを量子化テーブルQ1、CbとCrをQ2で量子化し、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;
}

struct YCbCr
{
	double Y, Cb, Cr;
};
struct RGB
{
	RGB(int r, int g, int b)
	{
		this->r = r;
		this->g = g;
		this->b = b;
	}
	int r, g, b;
};

YCbCr RGB2YCbCr(const RGB& rgb)
{
	YCbCr ycbcr;
	int R = rgb.r;
	int G = rgb.g;
	int B = rgb.b;

	ycbcr.Y = 0.299 * R + 0.5870 * G + 0.114 * B;
	ycbcr.Cb = -0.1687 * R - 0.3313 * G + 0.5 * B + 128;
	ycbcr.Cr = 0.5 * R - 0.4187 * G - 0.0813 * B + 128;

	return ycbcr;
}
RGB YCbCr2RGB(const YCbCr& ycbcr)
{
	double R, G, B;
	double Y = ycbcr.Y;
	double Cb = ycbcr.Cb;
	double Cr = ycbcr.Cr;
	R = Y + (Cr - 128) * 1.402;
	G = Y - (Cb - 128) * 0.3441 - (Cr - 128) * 0.7139;
	B = Y + (Cb - 128) * 1.7718;
	return RGB(R,G,B);
}

int main()
{
	PPM ppm("imori.pnm");
	int M = ppm.Get_width();
	int N = ppm.Get_height();
	std::vector<std::vector<double>> Y(M, std::vector<double>(N));
	std::vector<std::vector<double>> Cb(M, std::vector<double>(N));
	std::vector<std::vector<double>> Cr(M, std::vector<double>(N));
	for (int j = 0; j < N; j++)
		for (int i = 0; i < M; i++)
		{
			double r = ppm(i, j, 'r');
			double g = ppm(i, j, 'g');
			double b = ppm(i, j, 'b');
			YCbCr ycbcr = RGB2YCbCr(RGB(r, g, b));
			Y[i][j] = ycbcr.Y;
			Cb[i][j] = ycbcr.Cb;
			Cr[i][j] = ycbcr.Cr;
		}

	double q1[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 q2[64] = { 17, 18, 24, 47, 99, 99, 99, 99,
			   18, 21, 26, 66, 99, 99, 99, 99,
			   24, 26, 56, 99, 99, 99, 99, 99,
			   47, 66, 99, 99, 99, 99, 99, 99,
			   99, 99, 99, 99, 99, 99, 99, 99,
			   99, 99, 99, 99, 99, 99, 99, 99,
			   99, 99, 99, 99, 99, 99, 99, 99,
			   99, 99, 99, 99, 99, 99, 99, 99 };
	double Q1[8][8], Q2[8][8];
	for (int j = 0; j < 8; j++)
		for (int i = 0; i < 8; i++)
		{
			Q1[i][j] = q1[i + 8 * j];
			Q2[i][j] = q2[i + 8 * j];
		}

	PPM ppm2(M, 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>> iny(w, std::vector<double>(h));
			std::vector<std::vector<double>> incb(w, std::vector<double>(h));
			std::vector<std::vector<double>> incr(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++)
				{
					iny[i - w * I][j - h * J] = Y[i][j];
					incb[i - w * I][j - h * J] = Cb[i][j];
					incr[i - w * I][j - h * J] = Cr[i][j];
				}
			auto outy = dct2d(iny, w, h);
			auto outcb = dct2d(incb, w, h);
			auto outcr = dct2d(incr, w, h);
			for (int i = 0; i < w; i++)
				for (int j = 0; j < h; j++)
				{
					outy[i][j] = std::round(outy[i][j] / Q1[i][j]) * Q1[i][j];
					outcb[i][j] = std::round(outcb[i][j] / Q2[i][j]) * Q2[i][j];
					outcr[i][j] = std::round(outcr[i][j] / Q2[i][j]) * Q2[i][j];
				}
			auto iny2 = idct2d(outy, w, h, 8);
			auto incb2 = idct2d(outcb, w, h, 8);
			auto incr2 = idct2d(outcr, w, h, 8);
			for (int i = w * I; i < w * I + w; i++)
				for (int j = h * J; j < h * J + h; j++)
				{
					double y = iny2[i - w * I][j - h * J];
					double cb = incb2[i - w * I][j - h * J];
					double cr = incr2[i - w * I][j - h * J];
					YCbCr ycbcr;
					ycbcr.Y = y;
					ycbcr.Cb = cb;
					ycbcr.Cr = cr;
					RGB rgb = YCbCr2RGB(ycbcr);
					ppm2(i, j, 'r') = rgb.r;
					ppm2(i, j, 'g') = rgb.g;
					ppm2(i, j, 'b') = rgb.b;
					if (ppm2(i, j, 'r') < 0) ppm2(i, j, 'r') = 0;
					else if (ppm2(i, j, 'r') > 255) ppm2(i, j, 'r') = 255;
					if (ppm2(i, j, 'g') < 0) ppm2(i, j, 'g') = 0;
					else if (ppm2(i, j, 'g') > 255) ppm2(i, j, 'g') = 255;
					if (ppm2(i, j, 'b') < 0) ppm2(i, j, 'b') = 0;
					else if (ppm2(i, j, 'b') > 255) ppm2(i, j, 'b') = 255;
				}
		}
	ppm2.Flush("out.pnm");

	return 0;
}
```
![imori.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520857/4d9a9fb0-1d2a-0bc6-2fa5-fe1468914c45.png) ![out.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520857/0b1bf339-02de-2599-c9dd-1cfa8b130801.png)

フーリエ変換関連の問題がずいぶん続きましたが、これで一区切りです。