2007-09-13

DCT3回目開始！まずは、公式から

wikipedia等によると、公式はいくつかあるようなのですが、私は下記の公式を使う事にしました。
これから、以下の公式通りにプログラムを組んでみて、変換⇒逆変換と動かしてみて、
この公式が本当に正しいのかどうかを確かめてみようと思います。

===== 公式 =========================================================

N:変換する対象の要素数を表します。
π:円周率のパイです(3.1415926535・・・・)
S(x):DCT変換する前の値です。要素数N個の1次元の配列です。
S(u):DCT変換後の値です。要素数N個の1次元の配列です。
この値はDCT係数と呼ばれ、各々の波が、どの程度含まれているのか？という事を示しているそうです。

///// 変換 /////////////////////////////////////////////////////////

                        N-1                            
S(u) = √2/√N * C(u) * Σ  S(x) * cos( (2x+1)u*π/2N)  
                        x=0                            
但し、C(u)は、
    u=0  の時、1/√2
    u≠0 の時、1

///// 逆変換 ///////////////////////////////////////////////////////

                 N-1                                  
S(x) = √2/√N * Σ  C(u) * S(u) * cos( (2x+1)u*π/2N) 
                 u=0                                  
但し、C(u)は、
    u=0  の時、1/√2
    u≠0 の時、1     です。

====================================================================

↑が公式です。

今回は、要素の数を8個として実験します。
つまり、高さ1、幅8 の画像を1次元DCTするというイメージで行きたいと思います。

@	0	1	2	3	4	5	6	7
	10	20	30	40	50	60	70	80

次に、要素数を８としたときの公式と、そのプログラムを書きます。
上の公式の N が 8 になっただけです。

 /*
======変換=========================================
                     7
S(u) = 1/2 * C(u) * Σ  * S(x) * cos( (2x+1)uπ/16 )
                    x=0
但し、C(u)は、
    u=0  の時、1/√2
    u≠0 の時、1     です。
===================================================
*/
 double Cu = 0.0;
 double sigma = 0.0;
 double sigmain = 0.0;

 double cos=0.0;
 double cosin = 0.0;
 sb->AppendFormat("-----変換-----\r\n");

 for(u=0; u<8; u++){
 sigma = 0;
	 for(x=0; x<8; x++){
		 cosin = (2*x+1)*u*Math::PI/16;
		 cos = Math::Cos( cosin );
		 sigmain = Sx[x] * cos;
		 sigma += sigmain;
		 sb->AppendFormat("Sx[{4}]={5,4}, cos({0,8:F4})={1,8:f4},\tsigmain={2,10:f4},\tsigma={3,10:f4}\r\n"
			,cosin,cos,sigmain,sigma,x,Sx[x]);
	 }
	 if(u==0){
		 Cu = Math::Sqrt(1.0/2.0);
	 }else{
		 Cu = 1.0;
	 }
	 Su[u] = (int)(sigma/2.0*Cu);
	 sb->AppendFormat("Su[{0}]={1}\r\n",u,Su[u]);
 }

/*
=====逆変換========================================
              7
S(x) = 1/2 * Σ  * C(u) * S(u) * cos( (2x+1)uπ/16)
             u=0
但し、C(u)は、
	u=0  の時、1/√2
	u≠0 の時、1     です。
===================================================
*/
 sb->AppendFormat("-----逆変換-----\r\n");

 for(x=0; x<8; x++){
	 sigma = 0;

	 for(u=0; u<8; u++){
		 if(u==0){
			 Cu = Math::Sqrt(1.0/2.0);
		 }else{
			 Cu = 1.0;
		 }
		 cosin = (2*x+1)*u*Math::PI/16;
		 cos = Math::Cos(cosin);
		 sigmain = Su[u]*Cu*cos;
		 sigma += sigmain;
		 sb->AppendFormat("Su[{4}]={5,4}, cos({0,8:F4})={1,8:f4},\tsigmain={2,10:f4},\tsigma={3,10:f4}\r\n"
			,cosin,cos,sigmain,sigma,u,Su[u]);
	 }
	 Sx2[x] = (int)sigma/2;	 
	 sb->AppendFormat("Sx2[{0}]={1}\r\n",x,Sx2[x]);
 }
			 
 textBox3->Text = sb->ToString();

※計算途中の値を出力したいので、計算は、数回に分けてあります。

これを実行すると、以下のように出力され、この公式が、おそらく正しいであろうと思われるところまで来ました。

※double型からint型にキャスト(小数点切捨て)している為に、誤差が出ています。
また、DCT⇒逆DCTのみなら、プログラム上の丸め処理を考慮に入れなければ、可逆変換(完全に元に戻せる)です。

入力値	1	10	100	200	300	900	10	1
DCT後	538	-231	-417	420	-183	-119	357	-283
逆DCT後	1	10	100	199	299	899	10	0

入力値	100	100	100	100	100	100	100	100
DCT後	282	0	0	0	0	0	0	0
逆DCT後	99	99	99	99	99	99	99	99

2007-09-13

DCT2回目方針

JPEGでは、画像データを8x8ブロック毎に処理するため、2次元のDCT(離散コサイン変換)を使用していますが、
これ以降は、1次元のDCTでやってみようと思います。

1次元でも、2次元でも、数値を波に分解するという事には変わりないだろうと思っているからです。

プログラム環境は、Microsoft Visual C++ 2005 Express Edition , .Net Framework 2.0 です。
普通のC言語でも良いですが、自分の勉強もかねて、こっちにします。
http://www.microsoft.com/japan/msdn/vstudio/express/visualc/

また、グラフを描くソフトも使用します。
MATLABというソフトが有名なようですが、高価なようなので、フリーのMaTXというソフトにしました。
The MaTX Home Page

2007-09-13

DCT1回目注意事項

このブログに書いてある事は、正しいかどうかは分かりません。
このブログは、私にとって未知なDCTというものを、何とか理解してみようと試行錯誤するブログです。

私が間違っている事を書いているのを発見した場合は、教えていただけると幸いです。修正いたします。

また、私の数学力は、高校卒業(文系:数2、数B までです)

2007-09-13

DCT4回目 DCT係数を計算する所

///// 変換 /////////////////////////////////////////////////////////

                      7
Su(u) = 1/2 * C(u) * Σ  * Sx(x) * cos( (2x+1)uπ/16 )
                     x=0
但し、C(u)は、
    u=0  の時、1/√2
    u≠0 の時、1     です。

まず、1番目の波(DC成分と呼ばれ、実際には波ではなく、直線です)の成分、Su[0]の値を出す所。
上記の公式では、u = 0 となります。

よって、

                      7
Su(0) = 1/2 * 1/√2 * Σ  * Sx(x) * cos( (2x+1)*0*π/16 )
                     x=0

となりますので、
Σの中のcos() の中はゼロとなり、結果的に、

Su[0] = 1/2( 1/√2*Sx[0] + 1/√2*Sx[1] + ・・・・・・ + 1/√2*Sx[7] )

となります。
この値は、1/√2 に対して各要素が、どういう位置関係に居るのか？という値を個別に出していって、足したものでしょう。
この値を、全部足してしまう所に、私の不思議があります。
各要素に対する値を個別に保存しておくのならば、元に戻せるというのも納得です。
ですが、実際には全部足してしまうので、平均値(のような値)しか残りません。

元々の値が、
例１）100,100,100,100,100,100,100,100 というように一定の値であれば、足してしまっても構いませんが、
例２）0,100,100,100,100,100,100,200, というような値であっても、例１）と計算結果は同じになります・・・

※だから、８つの波で同様に値を出すのでしょうけど・・・

まぁ、次へ進みます

次は、2番目の波へ進みます。
同様に、公式に u=1 を代入すると以下となります。

              7
Su(1) = 1/2 * Σ Sx(x) * cos( (2x+1)*π/16 )
             x=0

ここで、cos() の中の値は、以下の様に変化します。
π*1/16, π*3/16, π*5/16, π*7/16, π*9/16, π*11/16, π*13/16, π*15/16

これは、ラジアンという単位で表されているようなので(初耳です)、高校時代の記憶通りに、(度：°)に変換すると
11.25°, 33.75°, 56.25°, 78.75°, 101.25°, 123.75°, 146.25°, 168.75°

これは、360°を16分割した22.5°区切りの真ん中の値かな？
0°-22.5°,22.5°-45°,45°-77.5°,77.5°-90°,90°-112.5°,112.5°-135°,135°-157.5°,157.5°-180°
まぁ、これは置いておいて

cos()を実際の値に直してみると

0.9808, 0.8315, 0.5556, 0.1951, -0.1951, -0.5556, -0.8315, -0.9808

となります。

この値で一つ感じるのは、
0番目と7番目、1番目と6番目、2番目と5番目、3番目と4番目、が、正負が逆になっています
なので、値が一定なら打ち消しあって 0 となります。

この事から、入力値が一定なら、Su[0] = 何らかの値、Su[1]、Su[2]、・・・ = 0 となるので、
元の値が復元出来そうな事は、何となく感じました。
ただ、以前、入力値がバラバラの時も復元出来るという所は不明です。

まぁ結局、計算は

Su(1) = 1/2( 0.9808*Sx[0] + 0.8315*Sx[1] + ・・・・・・ + -0.8315*Sx[6] + -0.9808*Sx[7])

となります。
これも、やっている事は、Su(0)の時と変わらず、
0.9808, 0.8315, 0.5556,・・・・というのは、波を表しているので、
これに対して、対応する要素の値との位置関係を求めているだけでしょう。

これも、個別の値を保持しているのなら分かりますが・・・全部足してしまうので、
結局残る値は、全体的に見てどうか？という事ですよね？？？？

うーーーん。
これ以降の、弟３、弟４・・・の波に関しても、波の波形が変わるだけなので、以下に表とグラフを載せておきます。

【コサイン内の角度】

	x=0	x=1	x=2	x=3	x=4	x=5	x=6	x=7
u=0	0°	0°	0°	0°	0°	0°	0°	0°
u=1	11.25°	33.75°	56.25°	78.75°	101.25°	123.75°	146.25°	168.75°
u=2	22.5°	67.5°	112.5°	157.5°	202.5°	247.5°	292.5°	337.5°
u=3	33.75°	101.25°	168.75°	236.25°	303.75°	11.25°	78.75°	146.25°
u=4	45°	135°	225°	315°	45°	135°	225°	315°
u=5	56.25°	168.75°	281.25°	33.75°	146.25°	258.75°	11.25°	123.75°
u=6	67.5°	202.5°	337.5°	112.5°	247.5°	22.5°	157.5°	292.5°
u=7	78.75°	236.25°	33.75°	191.25°	348.75°	146.25°	303.75°	101.25°

【コサインの実数値】

	x=0	x=1	x=2	x=3	x=4	x=5	x=6	x=7
u=0	0.7071	0.7071	0.7071	0.7071	0.7071	0.7071	0.7071	0.7071
u=1	0.9808	0.8315	0.5556	0.1951	-0.1951	-0.5556	-0.8315	-0.9808
u=2	0.9239	0.3827	-0.3827	-0.9239	-0.9239	-0.3827	0.3827	0.9239
u=3	0.8315	-0.1951	-0.9808	-0.5556	0.5556	0.9808	0.1951	-0.8315
u=4	0.7071	-0.7071	-0.7071	0.7071	0.7071	-0.7071	-0.7071	0.7071
u=5	0.5556	-0.9808	0.1951	0.8315	-0.8315	-0.1951	0.9808	-0.5556
u=6	0.3827	-0.9239	0.9239	-0.3827	-0.3827	0.9239	-0.9239	0.3827
u=7	0.1951	-0.5556	0.8315	-0.9808	0.9808	-0.8315	0.5556	-0.1951

【コサインの値に対するグラフ(実際の値を繋げたグラフ)】
DCTは、入力された値を、このグラフの各波に分解していると思います。

【コサインの値に対するグラフ(滑らかなグラフ)】

【各波を1枚に重ねたグラフ(実際の値を繋げたグラフ)】

【各波を1枚に重ねたグラフ(滑らかなグラフ)】

2007-09-11

一応、正常にデコードできるようになりました。

更新が遅れていましたが、前回のバグは修正でき、正常にデコードができるようになりました。
とりあえず、正常動作させる事を第一に書いていた為、人にお見せ出来るようなソースではありませんので、現段階では公開はしませんが、興味がある方がありましたら、おっしゃっていただければお送りいたします。

次回からは、私の本来の目的である、DCTとは？何をしているのか？という所を、突っ込んで書いていきたいと思います。

2007-09-09

ヤバイ！斬鉄剣出た

元画像はコチラ

あれから、いくつかの画像を正常に変換できていたので、
これで完成か？と笑っていたら甘かった・・・

2007-09-09

偶然縦縞、修正完了

テスト画像

修正前

修正後

ご覧のとおり、修正されました。

- - - - -

【豆知識】

JPEGは、(方式にもよりますが)画像をRGB⇒YCbCrという形式に変換してから圧縮します。
これが、ちょっとしたミソになっていて、
RGB(赤、緑、青)という色の表し方では、人間の目にとっては、RもGもBも、等価に価値のあるものですが、
YCbCrという色の表し方では、人間のからすると、Yは重要ですが、CbとCrは、そんなに重要ではありません。

つまり、RGBでは、どれか一つでも値が変わると、人間の目にとって色が変わりやすく
YCbCrでは、CbとCrの値が変わっても、人間の目にとっては、色が変わりにくいという事を表しています。

これを利用して、JPEGの中の圧縮率の高い方式では、
Yは、元の画像サイズと同じだけの量のデータを保存するが、
CbとCrは、元の画像よりも少ない量のデータしか保存しません(1/4とか)

例えば、
画像を、2ピクセル × 2ピクセルのデータ毎に見た場合に、
Yは、2x2で4ピクセル分のデータを保存
CbとCrは、(例えば)左上の1ピクセル分のデータしか保存しない
というような感じです。

これを、ダウンサンプリングと言うと思います。

- - - - -

そんなこんなで、今回のバグは、
ダウンサンプリングされている成分のデコードは上手く動いていたが(テスト画像ではCbCr)
ダウンサンプリングされていない成分のデコードの時に、格納する場所を間違えていたので、
8x8ブロック単位で考えた時に、Y方向に、2で割り切れるブロックの値を、-1ブロックした場所に上書きしていた
ので、空白(つまり 0x00 つまり黒っぽく)なっていたという事でした。