ガウシアンフィルタの並列処理実装の性能比較（リベンジ）

こんにてぃーは。前回ポストした記事にて、ガウシアンフィルタをSSEとOpenMPで並列処理してみたところ、SSEの効果が思ったほど良くなかったり、OpenCVのガウシアンフィルタ GaussianBlur() に比べて小カーネルでの性能がよろしくありませんでした。今日はそのリベンジであります。
とりあえずまずは「敵を知る」ということで、OpenCV2.3.1のソースコード*1を読んで調査しました。すげぇ苦労した甲斐あっていくつか知見が得られました。まずGaussianBlur()では、SSEによるベクトル演算化はされてますが、OpenMPなどようなマルチコア演算は未使用な模様。すごいですねぇ（マルチコア演算なしであの性能を超えれない＝男がすたるということですね）。・・・さらに細かく見ていくと、縦方向にフィルタする際、自分の実装のように転置→横フィルタ→転置とせずに、そのまま縦方向に処理していることが分かりました（まさかのCPUキャッシュガン無視っていう）。またちゃんとアラインメントを意識して高速化を図っていました。これはオレの設計がユニークすぎたというか、ただのアンポンタンだったようです。
これらの事実をふまえて、縦・横それぞれに関数を作ってみましょう。まずは横方向フィルタ処理から。

void convSseH(int w,int h,float* src,float* dst,float sigma)
{
    const int r=static_cast<int>(3.0*sigma); //３σハズレは無視
    const int r4=r-r%4+((r%4)?4:0); //rより大きい4の倍数を得る（アラインメント用）

    //カーネルの生成
    float* kernel=new float[r+1+r];
    const float ca=static_cast<float>(1.0/(sqrt(2.0*CV_PI)*sigma));
    const float cb=static_cast<float>(1.0/(2*sigma*sigma));
    for(int u=-r;u<=+r;u++)
        kernel[r+u]=ca*exp(-cb*u*u);

    //横方向のフィルタ処理
    #pragma omp parallel for
    for(int j=0;j<h;j++)
    {
        __m128 c,x0,s0;
        float* p=src+j*w;
        float* q=dst+j*w;
        for(int i=r4;i<w-r4;i+=4) //アライメント区切り毎に処理
        {
            s0=_mm_setzero_ps();
            for(int u=-r;u<=+r;u++)
            {
                c=_mm_load_ss(&kernel[r+u]);
                c=_mm_shuffle_ps(c,c,0);
                x0=_mm_loadu_ps(p+i+u);
                x0=_mm_mul_ps(x0,c);
                s0=_mm_add_ps(s0,x0);
            }
            _mm_store_ps(q+i,s0);
        }
    }
    delete kernel;
}

前回と大きく違うのは、カーネル中の一つの値を4画素に一気掛けして、同時に4画素分の出力画素値を得る点ですね。この構成だと、出力時にアラインメントを揃えた形で代入でき処理が高速になります。カーネル値も一つずつ取り出されるので、カーネルを保持する配列のアラインメントは気にしなくてOKです。意外とシンプルで短く書けますね。
次に新参の縦方向フィルタ処理です。

void convSseV(int w,int h,float* src,float* dst,float sigma)
{
    const int r=static_cast<int>(3.0*sigma);

    //カーネルの生成（対称性を利用するので右半分だけ用意）
    float* kernel=new float[1+r];
    const float ca=static_cast<float>(1.0/(sqrt(2.0*CV_PI)*sigma));
    const float cb=static_cast<float>(1.0/(2*sigma*sigma));
    for(int v=0;v<=r;v++)
        kernel[v]=ca*exp(-cb*v*v);

    //縦方向のフィルタ処理
    #pragma omp parallel for
    for(int j=r;j<h-r;j++)
    {
        __m128 c,x0,x1,s0;
        float* p=src+j*w;
        float* q=dst+j*w;
        for(int i=0;i<w;i+=4) //アライメント区切り毎に処理
        {
            c=_mm_load_ss(&kernel[0]);
            c=_mm_shuffle_ps(c,c,0);
            s0=_mm_load_ps(p+i);
            s0=_mm_mul_ps(s0,c);
            for(int v=1;v<=r;v++)
            {
                c=_mm_load_ss(&kernel[v]);
                c=_mm_shuffle_ps(c,c,0);
                x0=_mm_load_ps(p-v*w+i); //上側
                x1=_mm_load_ps(p+v*w+i); //下側
                x0=_mm_add_ps(x0,x1); //「折紙を半分に折って上と下くっつける」みたいな
                x0=_mm_mul_ps(x0,c);
                s0=_mm_add_ps(s0,x0);
            }
            _mm_store_ps(q+i,s0);
        }
    }
    delete kernel;
}

こちらも基本方針は一つのカーネル値について一気に4つ乗算です。フィルタはあくまで縦方向ですが、4つ同時に扱われるのは横方向の4つ（メモリが連続してるほう）です。これにより画素値はアラインメント区切り毎に処理される形にできます。またガウスカーネルの対称性より、係数が同じ部分（上側と下側）を加算してから乗算することで効率化を図ります。
で、早速ですが検証しましょう。今回は前回の反省を踏まえ、高解像度(4096×3072)な8bitグレースケール画像を用意し、σ=1〜128までテストすることにします。その他スペックもろもろは前回と同じです。

で、σとその計算時間のグラフです。すげー勢いで時間が増えていくので縦軸も対数軸にしました（ちょっと見辛くなりましたが）。またマルチコア演算のほうは、やっても前回と傾向が変わらなかったのと、OpenCVで未使用なので不公平という理由から今回載せてません。・・・結果の検証ですが、OpenCVの速度にはかなり追従するようになりました。が、まだ少し負けてますね。もう一つ重要な結果として、大カネールのときには単純な畳み込み実装のほうが速いことが確認できます。・・・これらの原因は明らかにアレなんですが、これについては次回のネタとしましょう。
ではまた近いうちにポストします。

*1:汚すぎてマ・ヂ・デ・ホ・ネ・ガ・オ・レ・タ・ヨ☆