画像変換のDCTをプログラム組んでみる その3 LLMのIDCT変換

三回目はLLMでのIDCT変換です。

これはさすがにバタフライ演算の形で書くのがいやのなので直接変換のコードを見てください。

コードはこんな感じです。

void LLMIDCT(float x[8],const float y[8])
{
	float a0,a1,a2,a3,b0,b1,b2,b3; float z0,z1,z2,z3,z4; float r[8]; int i;
	for(i = 0;i < 8;i++){ r[i] = (float)(cos((double)i / 16.0 * M_PI) * M_SQRT2); }
	z0 = y[1] + y[7]; z1 = y[3] + y[5]; z2 = y[3] + y[7]; z3 = y[1] + y[5];
	z4 = (z0 + z1) * r[3];
	z0 = z0 * (-r[3] + r[7]);
	z1 = z1 * (-r[3] - r[1]);
	z2 = z2 * (-r[3] - r[5]) + z4;
	z3 = z3 * (-r[3] + r[5]) + z4;
	b3 = y[7] * (-r[1] + r[3] + r[5] - r[7]) + z0 + z2;
	b2 = y[5] * ( r[1] + r[3] - r[5] + r[7]) + z1 + z3;
	b1 = y[3] * ( r[1] + r[3] + r[5] - r[7]) + z1 + z2;
	b0 = y[1] * ( r[1] + r[3] - r[5] - r[7]) + z0 + z3;
	z4 = (y[2] + y[6]) * r[6];
	z0 = y[0] + y[4]; z1 = y[0] - y[4];
	z2 = z4 - y[6] * (r[2] + r[6]);
	z3 = z4 + y[2] * (r[2] - r[6]);
	a0 = z0 + z3; a3 = z0 - z3;
	a1 = z1 + z2; a2 = z1 - z2;
	x[0] = a0 + b0; x[7] = a0 - b0;
	x[1] = a1 + b1; x[6] = a1 - b1;
	x[2] = a2 + b2; x[5] = a2 - b2;
	x[3] = a3 + b3; x[4] = a3 - b3;
	for(i = 0;i < 8;i++){ x[i] *= 0.5f; }
}

これがまた難しい形をしているのでわかりづらいですが、実際に計算してみれば逆変換となることがわかると思います。

プログラムコードとアルゴリズムについて要点を。

やっぱり係数は√2倍しておく

LLMの場合は正変換・逆変換ともに係数は√2倍したものを使います。おそらく実際に計算してみればどういうことかわかると思います。

バタフライ演算を行うとこの√2がきれいに消えるのがとてもうまいな~と思います。

逆変換時は一部係数の加減算が存在する

この部分が逆変換時の高速化の鍵で、この係数の計算をあらかじめ行っておけば高速になるのはすぐに理解できると思います。

もちろん、直値で組み込んでもかまいません。が、直値だとコードを見たときにどの係数の加減算だったかおそらくわからなくなると思います。

現に私もこのアルゴリズムを記録してあるノートにはこの形のまま書いていますが、プログラムコードは計算済みのものを使っているので

見た目何の係数だったかわからないです。はい。

cosの係数の性質を使って崩すことはできなくはないですが、逆変換時は計算しておいた定数を使うのが正しいと思います。

この係数の意味を出すのに計算された定数から逆変換するのがとても大変でした・・・。

最後の定数乗算部は画像処理では一つにまとめて

正変換と同様に画像処理では定数乗算部を二回以上行うのでまとめて行うようにします。

というわけで簡単にLLM変換を見てきました。このコードはすべて浮動小数で処理されていますが、

実際に動画圧縮などで使用するときは整数として処理した方が処理が早くなる傾向があります。

(特に整数の場合16bitである間はx8の処理であることと積和演算が使えること)

そちらは自分で考えてみてください。LLMをSSE化するときも同様です。その場合はこのコードをほとんどそのままSSEにすればいいだけだったりしますが。

コメントを残す

メールアドレスが公開されることはありません。 が付いている欄は必須項目です

この記事のトラックバック用URL