ここひと月、というか これ を作ってからアフィン変換に Bresenham のアルゴリズムを使えないか考えている。
理論上アフィン変換はメモリのコピーだけで、完全にメモリネック。つまりプリフェッチを行う以外に最適化の方法は無い。そして 1 ピクセルずつアクセスしている以上プリフェッチも難しい。(取り出してしまったキャッシュラインは全て処理する、とか考えられるけれども超複雑)
どころがどっこい素直に書いたアフィン変換の場合、ボトルネックはメモリではなく実は丸め処理にある。そこで丸めずに済む Bresenham の出番というわけ。
そんなわけで Bresenham の式を一般化しようとしてるんだけれども、中々うまい解が得られない。多分ググれば沢山あるんだろうけど、悔しいので探さない。
今のところ得られている式はこんな感じ。
x...イテレータ (int)
w...イテレートする範囲 (float)
h...イテレートする間に欲しい値の範囲 (float)
b...開始値 (float)
int(x) =floor(x)
round(x)=int(x+0.5)
基本的な一次方程式
y = round(ax+b)
= round(x*(h/w)+b)
// round(x)->int(x+1/2)
= int(x*(h/w)+b+1/2)
// 通分
= int(x*(2h/2w)+(2bw/2w)+(w/2w))
// 整数部と小数部を分離 (x*int(2h/2w))
= x*int(2h/2w) + int((2bw+w)/2w)
= x*int(h/w)+int(x*(2h%2w)) + int((2bw+w)/2w)
= x*int(h/w) + int({x*(2h%2w)+2bw+w}/2w)
// 整数部と小数部を分離 ((2bw+w)/2w)
= x*int(h/w)+int((2bw+w)/2w) + int({x*(2h%2w)+(2bw+w)%2w}/2w)
= x*int(h/w)+int((2b+1)/2) + int({x*(2h%2w)+(2bw+w)%2w}/2w)
// 最終形態
= x*int(h/w)+int((2b+1)/2)
+int({x*int(2h%2w)+int((2bw+w)%2w)}/int(2w))
ただこれ、イコールで結んじゃってるけれども元の式と同じ値にならない。1 だけ誤差が出る。
一つ分っているのは w や h が実数の時に int(2h%2w) とかやっちゃうと小数部分が死んでしまうところ。
精度を保つ必要があるのは 0 <= x < w の間だけなので、下のように分母に w を掛けてしまえば解決する。
y = x*int(h/w)+int((2b+1)/2)
+int({x*int(2hw%2ww)+int((2bw+w)w%2ww)}/int(2ww))
解決するんだけれども、w を二乗して更に 2 倍するので(入力値のビット数+1)*2ビットものビット数を用意しなくてはいけない。つまり 32 bit の CPU で処理する場合、入力値は 15 bit 以下である必要がある。
一般化したプログラムは以下の通り。
//! y = x*int(h/w)+int((2b+1)/2)+int({x*int(2hw%2ww)+int((2bw+w)w%2ww)}/int(2ww))
template<typename T=double, typename INT=int>
struct linear_iterator_round
{
INT y, dy, ydir;
INT f, df, flim;
linear_iterator_round( T w, T h, T start_y=0 )
: y(floor((2.0*start_y+1)*0.5)) // int((2b+1)/2)
, dy(floor(h/w)) // int(h/w)
, ydir(w*h != 0)
, f(fabs(floor(mod((2.0*start_y*w+w)*w, 2.0*w*w)))) // int((2bw+w)w%2ww)
, df(fabs(floor(mod(2.0*h*w, 2.0*w*w)))) // int(2hw%2ww)
, flim(fabs(floor(2.0*w*w))) // int(2ww)
{}
static T mod(T a, T b)
{
T r = fmod(a, b);
if ( a<0 != b<0 ) r += b;
return r;
}
//! 現在の値を得る
INT get() const { return y; }
//! 次の値を計算する
INT next()
{
y += dy;
f += df;
if ( f >= flim )
{
y += ydir;
f -= flim;
}
return y;
}
};
これを使うと floor(x+0.5) や round(x) (C99) を連発する場合に比べて 10 倍以上速くなる。
ただし最近の CPU は浮動小数点数の演算が速い為、これを利用する必要は無い。
C99 の lrint が使える場合、またはコンパイラの組み込み関数で浮動小数点->整数の変換が出来る場合(cvtxxx や __controlfp+intへのキャスト)はそちらを利用した方が速い。
floor がなんで遅いかって言うと、多分浮動小数点数から浮動小数点数に変換するから。浮動小数点数を整数レジスタに変換する場合は 1e30 とかって値を正しく変換出来ないけど、floor は変換後も double だから正しく元の値を保持しなくてはいけない。そして x86 や PPC なんかでは丸めた結果を浮動小数点数レジスタへ出力する命令は無いのでソフトウェア処理する羽目になるって寸法。
追記。ハードウェアで丸める場合も、古いタイプだと丸めモードをいちいち変更しなきゃいけないのが大きい。最近の x86 は丸めモード無しに直接丸める命令がある。
実用上 32,768 ピクセルを超える事も無いっちゃ無いんだけれども、もうちょっとひねって演算に必要なビット数を節約する方法は無いものかなあ。