Perlで離散コサイン変換

Perl

唐突にPerlのみでJpegを解析できるようなものを作ろうと思ったわけですが、挫折しそうです。
Jpegでは量子化のために離散コサイン変換（DCT）を行います。
通常の公式をそのまま使うと結構な時間がかかるので最適化を考えます。


参考文献

画像圧縮アルゴリズム - (6) JPEG法(1)
http://www2.starcat.ne.jp/~fussy/algo/algo8-6.htm


簡単な最適化の考え方としては、計算回数を減らすのと、同じ計算をしないことが挙げられます。


計算回数を減らすというのは、たとえば、
(x+1)(y-1) = (xy-1)
という簡単な因数分解（因数展開）ができる式を使うことになり、xとyの値が決まったとき、
左は加算1回、減算1回、乗算1回を行います。
右は乗算1回、減算1回となり、左とくらべ、加算1回分だけ短縮できたことになります。


同じ計算をしないこととは、単純に以前の計算結果を取っておくだけです。
たとえば、円周率PIを一々求めるのは面倒ですから、3.1415...という定数にしている場合が多いです。
今回のDCTの場合、cos部分の計算結果はN=8の時、8x8=64通りありますが、
計算回数を減らした場合でも、cos部分の計算は1024回行われます。
もし、この64通りをキャッシュできれば、960回は同じ計算をしなくて済む事になります。


またプログラミング言語が処理するのに得意な書き方があります。
単純なインクリメントでも、
$i=$i+1;
$i++;
の2種類ありますが、Perlでは後者のほうが早かったりします。今回では、

($dct_cos_table[$i][$y] ||= cos((2*$i+1)*$y*$T));

という辺りでそれを生かしています。同じ動作をする、

($dct_cos_table[$i][$y] ? $dct_cos_table[$i][$y] :($dct_cos_table[$i][$y] = cos((2*$i+1)*$y*$T)));

と比べると、コード量も少なく、実効速度も10%ほど早くなります。


use strict;
use warnings;
use Benchmark;

# 定数
my $M_PI = atan2(1,1)*4;
my $SQRT_HALF = 1/sqrt(2);

my @MCU = (
	[126, 138, 135, 118, 118, 126, 126, 130],
	[150, 168, 161, 122, 105, 109, 100, 118],
	[150, 150, 126, 150, 142, 126, 126, 117],
	[150, 161, 168, 130, 134, 150, 138, 130],
	[130, 118, 134, 142, 157, 142, 117, 126],
	[115, 117, 108, 117, 101,  99, 117, 126],
	[122, 130, 130, 138, 117, 108, 108, 138],
	[142, 118, 134, 117, 109,  91, 126, 109],
);
Benchmark::cmpthese(
	Benchmark::timethese(
		1000,
		{
		dct_idct_1 => sub {my (@dct,@idct); dct_1(\@MCU, \@dct);  idct_1(\@dct, \@idct);},
		dct_idct_2 => sub {my (@dct,@idct); dct_2(\@MCU, \@dct);  idct_2(\@dct, \@idct);},
		dct_idct_3 => sub {my (@dct,@idct); dct_3(\@MCU, \@dct);  idct_3(\@dct, \@idct);},
		}
	)
);

sleep 5;


sub dct_3{
	my ($ref_target, $ref_result) = @_;
	my $n = @$ref_target; 
	my $T = $M_PI/($n*2);
	my ($x,$y,$i);
	my @tmp=();
	my @dct_cos_table = ();

	for($x=0; $x<$n; ++$x){
		$ref_result->[$x] = [];
		$tmp[$x]=[];
		for($y=0; $y<$n; ++$y){
			$tmp[$x][$y]=0;
			
			for($i=0; $i<$n; ++$i){
				$tmp[$x][$y] += $ref_target->[$i][$y] *
				($dct_cos_table[$i][$x] ||= cos((2*$i+1)*$x*$T));
			}
			
			if($x==0){ $tmp[$x][$y] *= $SQRT_HALF; }
		}
		
		for($y=0; $y<$n; ++$y){
			$ref_result->[$x][$y] = 0;
			for($i=0; $i<$n; ++$i){
				$ref_result->[$x][$y] += $tmp[$x][$i] *
				($dct_cos_table[$i][$y] ||= cos((2*$i+1)*$y*$T));
			}
			
			if($y==0){ $ref_result->[$x][$y] *= $SQRT_HALF; }
			$ref_result->[$x][$y] *= (2/$n);
		}
	}
}

sub idct_3{
	my ($ref_target, $ref_result) = @_;
	my $n = @$ref_target;
	my ($x,$y,$i,$j);
	my $T = $M_PI/($n*2);
	my @dct_cos_table = ();
	my @tmp=();
	
	for($x=0; $x<$n; ++$x){
		$ref_result->[$x] = [];
		$tmp[$x]=[];
		for($y=0; $y<$n; ++$y){
			$tmp[$x][$y]=0;
			for($i=0; $i<$n; ++$i){
				$tmp[$x][$y] += ($i==0 ? $SQRT_HALF : 1) *
					$ref_target->[$i][$y] *
					($dct_cos_table[$x][$i] ||= cos((2*$x+1)*$i*$T));
			}
		}
		
		for($y=0; $y<$n; ++$y){
			$ref_result->[$x][$y] = 0;
			for($i=0; $i<$n; ++$i){
				$ref_result->[$x][$y] += ($i==0 ? $SQRT_HALF : 1) *
					$tmp[$x][$i] *
					($dct_cos_table[$y][$i] ||= cos((2*$y+1)*$i*$T));
			}
			$ref_result->[$x][$y] *= (2/$n);
		}
	}
}



sub dct_2{
	my ($ref_target, $ref_result) = @_;
	my $n = @$ref_target; 
	my ($x,$y,$i);
	my @tmp=();

	for($x=0; $x<$n; ++$x){
		$ref_result->[$x] = [];
		$tmp[$x]=[];
		for($y=0; $y<$n; ++$y){
			$tmp[$x][$y]=0;
			
			for($i=0; $i<$n; ++$i){
				$tmp[$x][$y] += $ref_target->[$i][$y] *
				cos((2*$i+1)*$x*$M_PI/(2*$n));
			}
			
			if($x==0){ $tmp[$x][$y] *= $SQRT_HALF; }
		}
		
		for($y=0; $y<$n; ++$y){
			$ref_result->[$x][$y] = 0;
			for($i=0; $i<$n; ++$i){
				$ref_result->[$x][$y] += $tmp[$x][$i] *
				cos((2*$i+1)*$y*$M_PI/(2*$n));
			}
			
			if($y==0){ $ref_result->[$x][$y] *= $SQRT_HALF; }
			$ref_result->[$x][$y] *= (2/$n);
		}
	}
}

sub idct_2{
	my ($ref_target, $ref_result) = @_;
	my $n = @$ref_target;
	my ($x,$y,$i,$j);
	my @tmp=();
	
	for($x=0; $x<$n; ++$x){
		$ref_result->[$x] = [];
		$tmp[$x]=[];
		for($y=0; $y<$n; ++$y){
			$tmp[$x][$y]=0;
			for($i=0; $i<$n; ++$i){
				$tmp[$x][$y] += ($i==0 ? $SQRT_HALF : 1) *
					$ref_target->[$i][$y] *
					cos((2*$x+1)*$i*$M_PI/(2*$n));
			}
		}
		
		for($y=0; $y<$n; ++$y){
			$ref_result->[$x][$y] = 0;
			for($i=0; $i<$n; ++$i){
				$ref_result->[$x][$y] += ($i==0 ? $SQRT_HALF : 1) *
					$tmp[$x][$i] *
					cos((2*$y+1)*$i*$M_PI/(2*$n));
			}
			$ref_result->[$x][$y] *= (2/$n);
		}
	}
}




sub dct_1{
	my ($ref_target, $ref_result) = @_;
	my $n = @$ref_target; # 幅
	
	for(my $x=0; $x<$n; ++$x){
		$ref_result->[$x] = [];
		for(my $y=0; $y<$n; ++$y){
			$ref_result->[$x][$y] = 0;
			for(my $i=0; $i<$n; ++$i){
				for(my $j=0; $j<$n; ++$j){
					$ref_result->[$x][$y] +=
						$ref_target->[$i][$j] *
						cos((2*$i+1)*$x*$M_PI/(2*$n)) *
						cos((2*$j+1)*$y*$M_PI/(2*$n));
				}
			}
			$ref_result->[$x][$y] *= (2/$n)*($x==0 ? $SQRT_HALF : 1)*($y==0 ? $SQRT_HALF : 1); # FA
		}
	}
	
}

sub idct_1{
	my ($ref_target, $ref_result) = @_;
	my $n = @$ref_target;
	
	for(my $x=0; $x<$n; ++$x){
		for(my $y=0; $y<$n; ++$y){
			$ref_result->[$x][$y] = 0;
			for(my $i=0; $i<$n; ++$i){
				for(my $j=0; $j<$n; ++$j){
					$ref_result->[$x][$y] +=
						($i==0 ? $SQRT_HALF : 1)*($j==0 ? $SQRT_HALF : 1) * 
						$ref_target->[$i][$j] *
						cos((2*$x+1)*$i*$M_PI/(2*$n)) *
						cos((2*$y+1)*$j*$M_PI/(2*$n));
				}
			}
			$ref_result->[$x][$y]*=(2/$n); # FA
		}
	}
}


dct_idct_1: 15 wallclock secs (14.48 usr +  0.00 sys = 14.48 CPU) @ 69.04/s (n=1000)
dct_idct_2:  3 wallclock secs ( 2.56 usr +  0.00 sys =  2.56 CPU) @ 390.32/s (n=1000)
dct_idct_3:  2 wallclock secs ( 1.84 usr +  0.00 sys =  1.84 CPU) @ 542.30/s (n=1000)
             Rate dct_idct_1 dct_idct_2 dct_idct_3
dct_idct_1 69.0/s         --       -82%       -87%
dct_idct_2  390/s       465%         --       -28%
dct_idct_3  542/s       685%        39%         --


DCTをしてiDCTをする動作を1000回繰り返して時間をみる。
結果的に9倍まで早くなった。


dct_idct_1は変換の公式をそのままコードに落とした感じで、
dct_idct_2はcos部分の計算回数が少なく済むようにした式を使ったもの（参考文献より）
dct_idct_3はdct_idct_2に加えて、cosの計算結果を記憶しているだけ。


効率がよくなるほど、コード量が多くなりかつメモリも消費するのは仕方がないね。
C言語で書けば1000/s行くのだろうか。