MATLABを使った離散コサイン基底と離散コサイン変換入門

No.2 画像データ圧縮篇

Version 1.2
　
早稲田大学教育・総合科学学術院　新井仁之
　
　
概要　2次元データ（主に画像データ）の2D離散コサイン変換を使った解析，特に画像データ圧縮やJPEG画像データ圧縮の離散コサイン変換を使った部分を解説する．
clear

１．画像データと線形空間

はじめに画像データを扱う上で必要な線形代数の基本的概念を説明しておく．
ピクセル（画素）のディジタル画像は，数値の M 行 N 列配列である．すなわち
(A1)
と表される配列である．ディジタル画像の場合，各成分 は の整数であるが，諸操作を加えると，一般には実数，場合によっては複素数になる．本稿では，実数からなる配列 (A1) 全体のなす集合を と表すことにする．
また (1) の配列を記述を簡約するため
と表す．特にすべての成分が 0 である 配列を または と表し，零元という．　
に対しては次の線形演算と内積を定義しておく．
線形演算
に対して
(和)
(スカラー積)
と定義する．
内積
に対して
とする．特に であるとき， と は直交するという．
次は内積のプログラム：
function C = InnerP2(A,B)
% A,B は共に (p,q)-行列．
C = sum(A.*B,'all');
end
基底と正規直交基底
基底と正規直交基底も の場合と同様に次のように定義される．なお の正規直交基底については [2] を参照．
定義1 （基底）
が の一組の基底であるとは，次の (B1), (B2) が成り立つことである．
(B1) 任意の は，ある実数 を選んで，必ず
　　　
　　と表すことができる．
(B2) は線形独立である．すなわち， をみたす実数 は 0 に限る．
コメント　なお， が の基底ならば， になることが証明される．すなわち の基底をなすベクトルの個数は MN 個である必要がある．これは が MN 次元の線形空間であることと関連している（次元の定義は[3]参照）．
定義2
が の一組の基底であり，かつ
(O1)
をみたすとき， の一組の正規直交基底であるという．
次の定理は正規直交基底であるかどうかを判定する際に便利である．
定理1
が (O1) をみたすならば， は の一組の正規直交基底である．

2. ベクトルのテンソル積と画像

二つの縦ベクトルから配列 (A1) を作る操作にテンソル積がある．
縦ベクトル と に対して， を の転置として
=
によりテンソル積 を定義する．
次のことが知られている（証明は [1] 参照）．
定理2
を の正規直交基底（定義は [2] 参照）とし， を の正規直交基底とする．このとき， に対して は の正規直交基底になっている．
この定理から， が の離散コサイン基底（定義は [2] に記してあるが，次節で再度述べる）であり， を の離散コサイン基底ならば， は の正規直交基底である．これを の 2D離散コサイン基底という．

3. 2D離散コサイン基底

まず の離散コサイン基底（こちらは1D離散コサイン基底という）の定義を記しておく．
,
に対して
.
の離散コサイン基底の定義も同様であるが，記しておこう．
,
に対して
.
これらのテンソル積として 2D離散コサイン基底が得られる．すなわち
定義3
に対して
を の2D離散コサイン基底という．
記号の簡略化のため，本稿では
と表すことにする．
離散コサイン基底を生成するプログラム
まず二つのベクトルのテンソル積を作る．
function C = TensorP(A,B)
%
% This function gives us the tensor product C of the following matrices A and B.
% A is (p,p)-matrix such that {A(:,1),...,A(:,p)} is a basis of R^p.
% B is (q,q)-matrix such that {B(:,1),...,B(:,q)} is a basis of R^q.
% C is (p,q,p,q)-matrix such that {C(:,:,j,k): j=1:p,k=1:q} is a basis of R^{pxq}.
%
p=size(A,1);
q=size(B,1);
C=zeros(p,q,p,q);
for j=1:p
for k=1:q
C(:,:,j,k)=A(:,j)*B(:,k)';
end
end
end
1D離散コサイン基底のテンソル積で2D離散コサイン基底を作る．
まずは [2] で作った1D離散コサイン基底が元になる：
function V = DCB(N)
% N 次元離散コサイン基底 V を行列形式で生成する．すなわち
% V は NxN 行列で，V(:,1),...,V(:,N) は R^N の離散コサイン基底．
V = zeros(N,N);
V(:,1)=1/sqrt(N)*ones(N,1);
for k = 2:N
for n = 1:N
V(n,k) = sqrt(2/N)*cos(pi/(2*N)*(2*(n-1)+1)*(k-1));
end
end
end
次に と の離散コサイン基底を作り，そのテンソル積で2D離散コサイン基底が作る．
TensorDCT(M,N) で の 2D離散コサイン基底が MxNxMxN 行列の形式で得られる（この方法だとメモリをかなり必要とするので注意）．たとえば v=TessorDCT(M,N) とすると，j=1,...,M, k=1,..,N に対して v(:,:,j,k) が の2D離散コサイン基底 となっている．
function DCB2 = TensorDCB(M,N)
A=DCB(M); %R^M の1D離散コサイン基底を表す MxM 行列
B=DCB(N); %R^N の1D離散コサイン基底を表す NxN 行列
DCB2 = TensorP(A,B);
end
例として の2D離散コサイン基底を作り，その強度画像を示す．これにより2D離散コサイン基底の様子がある程度わかる．以下では， を表している．
N = 8;
DCB88 = TensorDCB(N,N);
for j=1:N
for k=1:N
subplot(N,N,N*(j-1)+k)
imshow(DCB88(:,:,j,k),[])
title(['{\bf v}(',num2str(j-1),',',num2str(k-1),')'])
end
end
% saveas(gca, 'output.jpg') %画像を jpg 形式で保存したい場合．
2D離散コサイン基底の意味をもう少し詳しく述べておく．画像処理では「空間周波数」という考え方がある．詳細は説明しないが，要するに縞模様（あるいは格子模様）が細かいときに「高周波」，荒いときに「低周波」という．特に変化がなく一定の濃淡のものは「直流成分」という．
DCB88では左上が直流成分で，右下に行くほど低周波から高周波になっていく．この特性はデータ圧縮やノイズ低減で重要な役割を果たす．詳細は後述する．
Figure.jpg
図1
DCB88の正規直交性テスト
ここで得た の離散コサイン基底の正規直交性をMATLABを使って検証する．（あくまでも数値計算的に．厳密な証明は，たとえば[1]を参照）．
以下はDCT88の場合である．
%正規直交テスト行列 DeltaI とテスト行列 TestO の誤差が十分小なら数値的にOK.
M=8;N=8; %DCT88の検証，
DCBMN = TensorDCB(M,N);
TestO = zeros(M,N,M,N);
DeltaI = zeros(M,N,M,N);
for j=1:M
for k=1:N
for l=1:M
for m=1:N
TestO(j,k,l,m)=InnerP2(DCBMN(:,:,j,k),DCBMN(:,:,l,m));
if j==l && k==m
DeltaI(j,k,l,m)=1;
else
DeltaI(j,k,l,m)=0;
end
end
end
end
end
error = sum(abs(TestO - DeltaI),"all")
error = 4.6066e-13

4. 離散コサイン変換と逆離散コサイン変換

2D離散コサイン変換（2D DCT）と，2D逆離散コサイン変換（2D IDCT）を定義する準備ができた．
2D離散コサイン変換は， に対して，その周波数表現 を対応させる写像である．具体的には次の図の によって定義される．ただし下図では，記号の簡略化のため， と略記している．
2D逆離散コサイン変換は下図の である．
Decomp-Syntheesis_DCT_MATLAB3.jpg
この定義に基づき離散コサイン変換と逆離散コサイン変換をプログラムする．
[！] ただし，以下は理論的な考え方をプログラミングしたもので，高速アルゴリズムは考慮していない．
function DCT2 = DCTrans2(x)
% 2D離散コサイン変換（数学理論に沿ってプログラムしたもので，高速計算ではない）
[s1,s2] = size(x);
v = TensorDCB(s1,s2);
DCT2 = zeros(s1,s2);
for j = 1:s1
for k = 1:s2
w = v(:,:,j,k);
DCT2(j,k) = InnerP2(x,w);
end
end
end
function IDCT2 = IDCtrans2(x)
% 2D逆離散コサイン変換（数学理論に沿ってプログラムしたもので，高速計算ではない）
[s1,s2] = size(x);
v = TensorDCB(s1,s2);
IDCT2 = zeros(s1,s2);
for j = 1:s1
for k = 1:s2
IDCT2=IDCT2+x(j,k)*v(:,:,j,k);
end
end
end
簡単な例で検算，及びMATLABの dct2 との比較
x = magic(4) %テスト用データ
x = 4×4
16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1
まず本稿の数学理論重視に沿って作った（高速計算を考慮してない）プログラムの場合．
tic;Dx = DCTrans2(x);toc
経過時間は 0.001577 秒です。
MATLAB の　dct2
tic;MDx = dct2(x);toc
経過時間は 0.000539 秒です。
両者の誤差
sum(abs(Dx - MDx),'all')
ans = 4.3243e-14
反転公式の確認
IDCtrans2(Dx) %magic(4) が得られていることを確認．
ans = 4×4
16.0000 2.0000 3.0000 13.0000 5.0000 11.0000 10.0000 8.0000 9.0000 7.0000 6.0000 12.0000 4.0000 14.0000 15.0000 1.0000

5. 離散コサイン変換による画像データ圧縮事始め

離散コサイン変換を用いた画像データ圧縮の基本的な考え方は次のものである．
Figure2.jpg
ここで周波数の操作としては，単純なものとしては高周波をカットして周波数データ数を削減するものがある．まずはこれについて実験してみよう．簡単のため本稿ではグレースケールの画像を扱う．カラー画像については別途解説することにしたい．
データ圧縮の基本的な方法の一つは低周波は活かし，高周波は低減するというものである．
そこで，低周波の部分（上三角形の部分）は 1 で，高周波の部分（それ以外の部分）が 0 の配列を作成する．
function T = TriMatrix(M,N,R)
T = zeros(M,N);
if R>=N+M-1
T = ones(M,N);
else
for i = 1:M
for j = 1:N
if i+j<=R+1
T(i,j)=1;
end
end
end
end
end
TriMatrix(8,8,3) と TriMatrix(8,8,12) は次のようになる
TriMatrix(8,8,3)
ans = 8×8
1 1 1 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
TriMatrix(8,8,12)
ans = 8×8
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 0 1 1 1 1 1 0 0 0
Aを s1 x s2 画素の原画とする．T = TrigMatrix(s1,s2,R) とし，
とする．これが T により高周波をカットしたデータ圧縮画像になっている．
[！] なお，離散コサイン変換として，DCTrans2，IDCTrans2 を使うと，画像の画素数個の基底による展開を計算をすることになり，計算量や必要とするメモリーが膨大になる．
そのため以下ではMATLAB の dct2 と idct2 を使う
A_uint8 = imread('Photo_Gray.tif'); %本稿では筆者撮影の画像を使っているが，画像は好きなグレースケール画像を使う
A = double(A_uint8);
[s1,s2] = size(A)
s1 = 1200
s2 = 1600
G = [30 100 180]; % 1にする部分の行数
A_compress = zeros(s1,s2,length(G));
 
DCT_A = dct2(A);
 
for k = 1:length(G)
Q = TriMatrix(s1,s2,G(k));
A_compress(:,:,k) = idct2(Q.*DCT_A);
end
figure
subplot(2,2,1)
imshow(A,[])
title('Original Image')
for k = 1:length(G)
subplot(2,2,k+1)
imshow(A_compress(:,:,k),[])
title(['A_{',num2str(G(k)),'}'])
end
はほとんど概形しかわからないが， になると見た目の差はそれほど大きくない．本来ならば A の周波数表現の係数が 1200x1600=1920000 個必要であるにもかかわらず，0は勘定に入れなければ，では 180x180/2=16200 個
のみを使っている．

6. ブロック化した画像データ圧縮

JPEG形式では，離散コサイン変換を画像A全体に施すのではなく，画像を8x8の小さいブロックに分け，各ブロックに離散コサイン変換と周波数整形を施す方法を使っている．そこで，ここでは手始めにブロック化して TriMatrix を使った画像データ圧縮をしてみることにする．JPEG方式については次節で解説する．
MATLABには，画像をブロック化して，各ブロックに操作を施すコマンドに blockproc がある．これは
blockproc(画像，ブロックサイズ，ブロックに施す操作)
という構造になっている．いまブロックサイズは BlockSize = [8, 8] とし，これに dct2 を施すので，その操作を
fun = @(block_struct) dct2(block_struct.data)
とする．さらにこれを，例えば T = TriMatrix(8,8,1)で整形してから，逆離散コサイン変換することにすれば
fun = @(block_struct) idct2(T.*dct2(block_struct.data))
である．
画像としては既に使用している A （1200x1600画素）を使う．
BlockSize = [8,8];
T = TriMatrix(8,8,1);
fun = @(block_struct) idct2(T.*dct2(block_struct.data));
tic;
A1 = blockproc(A,BlockSize,fun);
toc
経過時間は 2.065703 秒です。
ブロックごとならば DCTrans2 でも，メモリー不足にはならずに計算ができる（dct2，idct2を使ったものよりかなり時間はかかる）．
fun = @(block_struct) IDCtrans2(T.*DCTrans2(block_struct.data));
tic;
A2 = blockproc(A,BlockSize,fun);
toc
経過時間は 8.734091 秒です。
原画像 A と上で得た圧縮画像を図示する．A1 と A2 は同じになる．
figure
subplot(2,2,1)
imshow(A,[])
title('原画像')
subplot(2,2,3)
imshow(A1,[])
title('A1')
subplot(2,2,4)
imshow(A2,[])
title('A2')
さて，画像表示が小さいとデータ圧縮による画像劣化がわかりにくいので，一部を取り出して比較する，データ圧縮画像では，各8x8ブロックでは平均値をあてがっていることがわかる．
A_part= A(1:64,1:64);
A1_part = A1(1:64,1:64);
% A2_part = A2(1:64,1:64);
figure
subplot(1,2,1)
imshow(A_part,[])
title('原画像の一部')
subplot(1,2,2)
imshow(A1_part,[])
title('A1の一部')

7. JPEG方式（ただしDCTを使う部分のみ）

JPEGでは TriMatrix のような素朴なものではなく，次のような量子化テーブルと丸め込みを組み合わせてた方法をとっている（たとえば[4] 参照）．以下では [4] を参考にした．
% 量子化テーブル
q=[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];
まず画像 A の各ブロックから 128 を引いて，0を中心のデータとする．B=A-128.
次にそれを離散コサイン変換する. freq_image = dct2(B).
量子化テーブルを n 倍して，その逆数を掛け，丸める．freq_image_surj = round(freq_image./(q*n)).
それを q*n 倍して128を加算して丸めるて，逆離散コサイン変換．synth = round(idct2(freq_image_surj*(q*n)+128)).
さらに逆離散コサイン変換して，（先ほど 128を引いたので，ここで）128を加える．
圧縮率をコントロールするパラメータとして quality が使われる．quality を大きくすると圧縮率が増加する．そのため画像の劣化も増加する．
% quality
n=10; % とりあえず 10 に設定しておく．
% ブロック処理 (cf. [4])
freq_image = @(X) dct2(X.data-128); % 画像-128 の周波数表現
freq_surg_image = @(X) round(X.data./(q*n)); % 周波数整形
synth = @(X) round(idct2((X.data).*(q*n))+128); % 再構成
 
A_freq = blockproc(A,[8,8],freq_image); %ブロック処理による周波数表現
A_freq_surg = blockproc(A_freq,[8,8],freq_surg_image); %ブロック処理による周波数整形
A_compress = blockproc(A_freq_surg,[8,8],synth); %ブロック処理による再構成
A_compress_uint8= uint8(A_compress); %uint8 への変換
disp(['全体の成分数=',num2str(numel(A_freq))])
全体の成分数=1920000
zero_component = numel(A_freq)-nnz(A_freq_surg); % nnz は非零の個数
disp(['周波数整形で 0 になった成分の個数=',num2str(zero_component)])
周波数整形で 0 になった成分の個数=1849257
ratio = (zero_component)/numel(A_freq);
disp(['圧縮率=',num2str(ratio)])
圧縮率=0.96315
原画像とデータ圧縮画像を比較する．
figure
montage({A_uint8,A_compress_uint8})
title(['左図：原画像，右図：圧縮率=',num2str(ratio),',','quality =', num2str(n),'の圧縮画像'])
saveas(gca,['jpegcompress',num2str(n),'.jpg']) %画像を jpg 形式で保存したい場合
練習問題
quality を変えて，データ圧縮画像を観察せよ．特に quality を上げると，圧縮率はよくなるが，ブロック状のひずみがでてくる．これはブロックひずみと呼ばれている．
余談：triuを使った方法との比較
tic;
A_triu = ones(10,10);
R = triu(A_triu,3);
RL = fliplr(R)
RL = 10×10
1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 0 0 0 0 1 1 1 1 1 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
toc
経過時間は 0.007252 秒です。
tic;
A_Tri=TriMatrix(10,10,7)
A_Tri = 10×10
1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 0 0 0 0 1 1 1 1 1 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
toc
経過時間は 0.004624 秒です。

参考文献

[1] 新井仁之，フーリエ解析学，朝倉書店，2003.
[2] 新井仁之，MATLABを使った離散コサイン基底と離散コサイン変換入門，No.1 1次元データ篇，http://www.araiweb.matrix.jp/MATLAB/DCT_Tutorial_using_MATLAB1
[3] 新井仁之，線形代数 基礎と応用，日本評論社，2006.
[4] A. McAndrew, A Computational Introduction to Digital Image Processing, 2nd ed., CRC Press, 2016.
履歴
Version 1.1 2024/07/03 公開
Version 1.2 2024/07/04 minor changes
Copyright © 2024 by Hitoshi Arai.