【離散コサイン変換】二次元DCTと対称性

# 背景
[【離散コサイン変換】DCT-IIの式をぼんやり理解する](https://qiita.com/jajagacchi/items/e7b427c71257860f95e6)
[【離散コサイン変換】DFTを使ってDCTしてみる](https://qiita.com/jajagacchi/items/2e3c40ce7795ade5b87c)
の続きです。

# 使用したライブラリ
[【画像処理100本ノック】独自の画像入出力クラスを作る](https://qiita.com/jajagacchi/items/834e8425b84f94df8a49)
[【画像処理100本ノックに挑戦】独自の離散フーリエ変換クラスを作る](https://qiita.com/jajagacchi/items/0b5232aa8ae99b9b8114)

# 二次元DCTと対称性
2次元の場合、離散フーリエ変換同様にある方向に一次元DCTした後に直交する方向に一次元DCTすればOKです。問題はフィルタ処理です。フィルタ処理は対称性を維持するように行う必要があります。

DCT-IIの場合はそれほど難しくありませんが、DCT-IIIはやっかいです。そこで、DCT-IIを例に考えてみます。一次元で考えます。$N$個のデータ$x_n$があるとします。DFTを用いて計算するために、$x_n$を用いて$4N$個の要素からなる$y_n$を定義します。

符号の変換関数$s(n)$を定義します。これは

```math
s(2n+1)=s(4N-2n-1)=1\;(0\leq n<N)
```

で与えられ、ほかの点はゼロです。

次に添え字の変換関数$t(n)$を用意します。

```math
t(2n+1)=t(4N-2n-1)=n\;(0\leq n<N)
```

で定義されます。ほかの点の値は任意です。

このとき

```math
y_{mn} = s_x(m)s_y(n)x_{t_x(m)t_y(n)}
```

としてDFTすればDCT-IIが得られます。対称性を維持するためには関数$t$によって同じ$x$に写像される場合は同じように変化させます。DCT-IIIの場合も$s,t$の定義が変わるだけで同様です。

まずはDCT-IIの場合のコードを示します。

```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;
}

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

	std::vector<std::vector<std::complex<double>>> f(M, std::vector<std::complex<double>>(N));
	for (int j = 0; j < N; j++)
		for (int i = 0; i < M; i++)
		{
			int r = ppm(i, j, 'r');
			int g = ppm(i, j, 'g');
			int b = ppm(i, j, 'b');
			int y = (std::round)(0.2126 * r + 0.7152 * g + 0.0722 * b);
			f[i][j] = y;
		}
	std::vector<double> sx(4 * M), sy(4 * N), tx(4 * M), ty(4 * N);
	for (int n = 0; n < M; n++)
	{
		sx[n] = tx[n] = 0;
	}
	for (int m = 0; m < N; m++)
	{
		sy[m] = ty[m] = 0;
	}

	for (int n = 0; n < M; n++)
	{
		sx[2 * n + 1] = sx[4 * M - 2 * n - 1] = 1;
		tx[2 * n + 1] = tx[4 * M - 2 * n - 1] = n;
	}
	for (int m = 0; m < N; m++)
	{
		sy[2 * m + 1] = sy[4 * N - 2 * m - 1] = 1;
		ty[2 * m + 1] = ty[4 * N - 2 * m - 1] = m;
	}

	std::vector<std::vector<std::complex<double>>> in(4 * M, std::vector<std::complex<double>>(4 * N));
	std::vector<std::vector<double>> in2(M, std::vector<double>(N));
	for (int n = 0; n < 4 * M; n++)
		for (int m = 0; m < 4 * N; m++)
		{
			in[n][m] = sx[n] * sy[m] * f[tx[n]][ty[m]];
		}
	for (int n = 0; n < M; n++)
		for (int m = 0; m < N; m++)
		{
			in2[n][m] = f[n][m].real();
		}
	DFT2D dft2d(4 * M, 4 * N);
	auto out = dft2d.forward(in);
	auto out2 = dct2d(in2, M, N);

	for (int n = 0; n < M; n++)
		for (int m = 0; m < N; m++)
		{
			std::cout << n << " " << m << " " << out[n][m].real()/4. << " " << out2[n][m] << std::endl;
		}
	return 0;
}
```

DCT-IIIでは次のようになります。

```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++)
		{
			double a, b;
			for (int m = 0; m < M; m++)
				for (int n = 0; n < N; 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;
				}
		}
	return ret;
}

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

	std::vector<std::vector<std::complex<double>>> f(M, std::vector<std::complex<double>>(N));
	for (int j = 0; j < N; j++)
		for (int i = 0; i < M; i++)
		{
			int r = ppm(i, j, 'r');
			int g = ppm(i, j, 'g');
			int b = ppm(i, j, 'b');
			int y = (std::round)(0.2126 * r + 0.7152 * g + 0.0722 * b);
			f[i][j] = y;
		}
	std::vector<double> sx(4 * M), sy(4 * N), tx(4 * M), ty(4 * N);
	for (int n = 0; n < M; n++)
	{
		sx[n] = tx[n] = 0;
	}
	for (int m = 0; m < N; m++)
	{
		sy[m] = ty[m] = 0;
	}
	for (int n = 0; n < N; n++)
	{
		sx[n] = 1;
		tx[n] = n;
	}
	for (int n = N + 1; n <= 2 * N; n++)
	{
		sx[n] = -1;
		tx[n] = 2 * N - n;
	}
	for (int n = 2 * N + 1; n < 3 * N; n++)
	{
		sx[n] = -1;
		tx[n] = n - 2 * N;
	}
	for (int n = 3 * N + 1; n < 4 * N; n++)
	{
		sx[n] = 1;
		tx[n] = 4 * N - n;
	}

	for (int m = 0; m < M; m++)
	{
		sy[m] = 1;
		ty[m] = m;
	}
	for (int m = M + 1; m <= 2 * M; m++)
	{
		sy[m] = -1;
		ty[m] = 2 * M - m;
	}
	for (int m = 2 * M + 1; m < 3 * M; m++)
	{
		sy[m] = -1;
		ty[m] = m - 2 * M;
	}
	for (int m = 3 * M + 1; m < 4 * M; m++)
	{
		sy[m] = 1;
		ty[m] = 4 * M - m;
	}

	std::vector<std::vector<std::complex<double>>> in(4 * M, std::vector<std::complex<double>>(4 * N));
	std::vector<std::vector<double>> in2(M, std::vector<double>(N));
	for (int n = 0; n < 4 * M; n++)
		for (int m = 0; m < 4 * N; m++)
		{
			in[n][m] = sx[n] * sy[m] * f[tx[n]][ty[m]];
		}
	for (int n = 0; n < M; n++)
		for (int m = 0; m < N; m++)
		{
			in2[n][m] = f[n][m].real();
		}
	DFT2D dft2d(4 * M, 4 * N);
	auto out = dft2d.forward(in);
	auto out2 = dct2d(in2, M, N);

	for (int n = 0; n < M; n++)
		for (int m = 0; m < N; m++)
		{
			std::cout << n << " " << m << " " << out[2*n+1][2*m+1].real()/16. << " " << out2[n][m] << std::endl;
		}
	return 0;
}
```
どちらも、定義通り計算したものと一致します。ただ、FFTを用いていないので無駄にメモリと時間を浪費します。$4\times4=16$倍無駄でしょうか。100本ノックにDFT版を使うのは無駄なので定義通り計算することにします。