離散コサイン変換(DCT)をPythonで実装した[Python]

ダブルミーニングでテストが終わりました．春休みに入ったので今年の目標の「数学を勉強する」を早速実行してみます．将来勉強するはめになりそうなものから学ぼうと思うので，まずは主軸変換から勉強します．

フーリエ変換はnumpyやscipyで実装されているので，されていないであろう離散コサイン変換(DCT)を今日は実装してみました……と思っていたらscipy.fftpackにdctの実装ありました（公式サイト)．こちらはWikipediaで紹介されている基底関数DCT-I，DCT-II，DCT-IIIが実装されているようです．

1.離散コサイン変換(DCT)の概要

離散コサイン変換(Discrete Cosine Transform)はWikipediaさんを見る限りではいろんな基底関数があるようですが，今手元にある本『画像処理とパターン認識入門』では離散コサイン変換の基底関数を次のように定義しています．Wikipediaで言うところのDCT-IIですね．
$$\phi_k[i] = \begin{cases} \frac{1}{\sqrt{N}} \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,(k=0) \\ \sqrt{\frac{2}{N}} \cos {\frac{(2i+1)k\pi}{2N}} (k=1,2,...,N-1) \end{cases}$$

この基底関数ですがこんな感じのグラフになります．N=64，横軸がiに相当し上からk=0,1,2,...,9の時の基底関数のグラフです．この定義だとほぼ普通の三角関数ですね．

離散化された信号をf_iとすると，離散コサイン変換，離散コサイン逆変換の式は以下のとおりです．
$$C[k] = \sum_{i=0}^{N-1}f_{i}\phi_k[i] \,\,(k=0,1,...,N-1)\\ f_{i}=\sum_{k=0}^{N-1}C[k]\phi_{k}[i] \,\,(i=0,1,...,N-1)$$
とりあえずこいつを実装します．

2.離散コサイン変換(1次元)の実装例

まぁこんな感じです．関数phiは基底関数です．本で紹介されていた数式以外にWikipediaにあったDCT-IVもコメントアウトしてありますが一応実装してあります．
ベンチマークとして乱数で作った配列データを試しに変換・逆変換してみます．

	#coding: utf-8
	import numpy as np
	import matplotlib.pyplot as plt
	"""
	参考:
	[1]『画像処理とパターン認識入門』酒井幸市 著
	[2] scipy.fftpack.dct http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.fftpack.dct.html
	"""
	class DCT:
	def __init__(self,N):
	self.N = N # データ数．
	# 1次元離散コサイン変換の変換行列を予め作っておく
	self.phi_1d = np.array([ self.phi(i) for i in range(self.N) ])
	def dct(self,data):
	""" 1次元離散コサイン変換を行う """
	return self.phi_1d.dot(data)
	def idct(self,c):
	""" 1次元離散コサイン逆変換を行う """
	return np.sum( self.phi_1d.T * c ,axis=1)
	def phi(self,k):
	""" 離散コサイン変換(DCT)の基底関数 """
	# DCT-II
	if k == 0:
	return np.ones(self.N)/np.sqrt(self.N)
	else:
	return np.sqrt(2.0/self.N)*np.cos((k*np.pi/(2*self.N))*(np.arange(self.N)*2+1))
	# DCT-IV(試しに実装してみた)
	#return np.sqrt(2.0/N)*np.cos((np.pi*(k+0.5)/self.N)*(np.arange(self.N)+0.5))
	if __name__=="__main__":
	N = 100 # データ数を100とします
	dct = DCT(N) # 離散コサイン変換を行うクラスを作成
	x = np.random.random_sample(N) # N個の乱数データを作成
	c = dct.dct(x) # 離散コサイン変換を実行
	y = dct.idct(c) # 離散コサイン逆変換を実行
	# 元のデータ(x)と復元したデータ(y)をグラフにしてみる
	plt.plot(x,label="original")
	plt.plot(y,label="restored")
	plt.legend()
	plt.title("original data and restored data")
	plt.show()

view raw dct_1dim.py hosted with ❤ by GitHub

実行結果

まぁ問題なくデータが復元されてますね．乱数なので高周波数成分があるので圧縮は無理ですが．

3.離散コサイン変換(2次元)の概要

さて手元の本のタイトルには「画像処理」の文字が書かれています．すなわち二次元データに対する解説も書かれています（というかむしろそれが本題）．なのでpythonで二次元データに対しても実装してみます．

本によれば二次元のDCT係数F_k,lは次のように表されるとのこと．また画像の濃度値f_i,jはDCT係数を用いてこのように逆変換できるとのこと．
$$F_{k,l}=\sum_{j=0}^{N-1}\sum_{i=0}^{N-1}f_{i,j}\phi_k[i]\phi_l[j] \\ f_{i,j}=\sum_{l=0}^{N-1}\sum_{k=0}^{N-1}F_{k,l}\phi_k[i]\phi_{l}[j]$$
ここでこの基底ベクトル（基底画像）がどんな形をしているか見てみます．

この画像，離散コサイン変換（DCT)で検索すると出てくるやつですね．これらの画像は10x10ピクセルの画像を100個並べたものですが，二次元離散コサイン変換ではこの100個の画像をうまいこと足し合わせることで様々な画像を表現するというのがテンションが上ります．（つまりこの画像の中にあんな画像，こんな画像まで，ありとあらゆる画像が含まれていると言っても過言ではない……？）

4.離散コサイン変換(2次元)の実装例

では実装です．scipyにも2次元の離散コサイン変換は実装されてないのでちょっとテンションが上ります．まぁ2次元のフーリエ変換は実装されてますが．

サンプルとして適当な画像を変換・逆変換してみます．

	#coding: utf-8
	import numpy as np
	import matplotlib.pyplot as plt
	"""
	参考:
	[1]『画像処理とパターン認識入門』酒井幸市 著
	[2] scipy.fftpack.dct http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.fftpack.dct.html
	"""
	class DCT:
	def __init__(self,N):
	self.N = N # データ数．
	# 1次元，2次元離散コサイン変換の基底ベクトルをあらかじめ作っておく
	self.phi_1d = np.array([ self.phi(i) for i in range(self.N) ])
	# Nが大きいとメモリリークを起こすので注意
	# MNISTの28x28程度なら問題ない
	self.phi_2d = np.zeros((N,N,N,N))
	for i in range(N):
	for j in range(N):
	phi_i,phi_j = np.meshgrid(self.phi_1d[i],self.phi_1d[j])
	self.phi_2d[i,j] = phi_i*phi_j
	def dct(self,data):
	""" 1次元離散コサイン変換を行う """
	return self.phi_1d.dot(data)
	def idct(self,c):
	""" 1次元離散コサイン逆変換を行う """
	return np.sum( self.phi_1d.T * c ,axis=1)
	def dct2(self,data):
	""" 2次元離散コサイン変換を行う """
	return np.sum(self.phi_2d.reshape(N*N,N*N)*data.reshape(N*N),axis=1).reshape(N,N)
	def idct2(self,c):
	""" 2次元離散コサイン逆変換を行う """
	return np.sum((c.reshape(N,N,1)*self.phi_2d.reshape(N,N,N*N)).reshape(N*N,N*N),axis=0).reshape(N,N)
	def phi(self,k):
	""" 離散コサイン変換(DCT)の基底関数 """
	# DCT-II
	if k == 0:
	return np.ones(self.N)/np.sqrt(self.N)
	else:
	return np.sqrt(2.0/self.N)*np.cos((k*np.pi/(2*self.N))*(np.arange(self.N)*2+1))
	# DCT-IV(試しに実装してみた)
	#return np.sqrt(2.0/N)*np.cos((np.pi*(k+0.5)/self.N)*(np.arange(self.N)+0.5))
	if __name__=="__main__":
	N = 10 # データの次元は10x10とする
	dct = DCT(10) # 離散コサイン変換を行うクラスを作成
	# サンプル画像を作る
	img = np.array([
	[0,0,0,0,0,0,0,0,0,0],
	[0,0,1,1,1,1,1,1,0,0],
	[0,1,0,0,0,0,0,0,1,0],
	[0,1,0,0,0,0,0,0,1,0],
	[0,0,0,0,0,0,1,1,0,0],
	[0,0,0,0,1,1,0,0,0,0],
	[0,0,0,1,0,0,0,0,0,0],
	[0,0,1,0,0,0,0,0,0,0],
	[0,1,1,1,1,1,1,1,1,0],
	[0,0,0,0,0,0,0,0,0,0]
	])
	c = dct.dct2(img) # 2次元離散コサイン変換
	y = dct.idct2(c) # 2次元離散コサイン逆変換
	# 元の画像と復元したものを表示
	plt.subplot(1,2,1)
	plt.imshow(img,cmap="Greys")
	plt.title("original")
	plt.xticks([])
	plt.subplot(1,2,2)
	plt.imshow(y,cmap="Greys")
	plt.title("restored")
	plt.xticks([])
	plt.show()

view raw dct_2dim.py hosted with ❤ by GitHub

実行結果です．

画像がちゃんと元の状態に戻っているのがわかります．ちなみにDCT係数（周波数成分）はこんな感じでした．

概ね右下に行くほど色が薄くなっています．ここから低周波成分(左上）が多くを占めていることがわかります．
高周波成分を消していったものがこちらです．

DCT係数cについてc[i:,i:] = 0とした時に復元された画像です．i=3でだいたい認識できそうですね．結構データを圧縮できてるのではないでしょうか．まぁまぁうまく行って良かったです．