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

こんにてぃは。今日は、画像処理における超基本ツールの一つと言えるガウシアンフィルタについて、C++でSSE（SIMD演算）とOpenMP（マルチコア演算）を使ったときの性能を比較してみました。

（【2012/04/23追記】コード内に無駄なキャストが残っていたので除去しました。）
まず基本方針ですが、ガウスカーネルの幅は±3σまでとし、単精度(float)で処理します。またOpenMPを使って各ラインを独立にマルチコア処理させます。ガチで効率的にやるならIIRフィルタのほうがベターかと思いますが、今回は畳み込み演算そのものの高速化に焦点を当てたいので、そういうのは今回は扱ってません。
何の工夫もない畳み込み演算による実装は次のようになります。念のため対称性を利用して乗算を半分程度に省いてます。

void convH(int w,int h,float* src,float* dst,float sigma)
{
    //ガウスカーネルの生成
    const int r=static_cast<int>(3.0*sigma)
    float* kernel=new float[r+1];
    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 x=0;x<=r;x++)
        kernel[x]=ca*exp(-cb*x*x);

    //畳み込み演算(OpenMP)
    #pragma omp parallel for
    for(int j=0;j<h;j++)
    {
        float* p=src+j*w;
        float* q=dst+j*w;
        for(int i=r;i<w-r;i++)
        {
            float sum=kernel[0]*p[i];
            for(int t=1;t<=r;t++)
                sum+=kernel[t]*(p[i-t]+p[i+t]);
            q[i]=sum;
        }
    }
    delete kernel;
}

SSEも併用した実装は以下のとおり。こちらはループの最後の処理を簡単にするため、カーネル幅が4の倍数になるように詰め物をしてます。内積演算したい場合、SSE4.1にDPPSという専用命令がありますが全く速度差がなかったので、今回はより対象マシンの広いSSE1で実装しました。またSSE使用時はメモリのアラインメントが重要になりますが、今回の実験では簡単のために無視しています（カーネルのアラインメントに関する処理はコメントアウトしています・・・が、これも微々たる差だったという）。

void convSseH(int w,int h,float* src,float* dst,float sigma)
{
    //ガウスカーネルの生成
    const int r=static_cast<int>(3.0*sigma);
    const int sz=calcAlignedSize<4>(2*r+1);
    //float* kernel=reinterpret_cast<float*>(_aligned_malloc(sizeof(float)*sz,16));
    float* kernel=new float[sz];
    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 x=-r;x<=+r;x++)
        kernel[x+r]=ca*exp(-cb*x*x);
    for(int i=2*r+1;i<sz;i++)
        kernel[i]=0.0;

    //畳み込み演算(OpenMP+SSE)
    #pragma omp parallel for
    for(int j=0;j<h;j++)
    {
        __m128 regK,regA,sums;
        float* p=src+j*w;
        float* q=dst+j*w;
        for(int i=0;i<w-sz;i++)
        {
            sums=_mm_setzero_ps();
            for(int t=0;t<sz;t+=4)
            {
                regK=_mm_loadu_ps(&kernel[t]); //regK=_mm_load_ps(&kernel[t]);
                regA=_mm_loadu_ps(&p[i+t]);
                regA=_mm_mul_ps(regA,regK); //各成分を乗算
                sums=_mm_add_ps(sums,regA); //その乗算結果を累積
            }
            _declspec(align(16)) float buf[4]={0};
            _mm_store_ps(buf,sums); //累積結果をbufに転送
            q[i+r]=buf[0]+buf[1]+buf[2]+buf[3]; //各累積結果を一つにまとめる
        }
    }
    delete kernel;
    //_aligned_free(kernel);
}

・・・上記のコードは横方向のみのフィルタです。縦方向にもフィルタする場合は、一度データを転置してからもう一度横方向にフィルタリングし、その結果を再度転置します。メモリアクセスは局所的に連続したほうが高速なため、わざわざでも転置したほうが全体的には速いことが多いです。うまいことバッファを設ければ転置せずにさらに高速に処理できる可能性もありますが、面倒なのでそっち方面は考えないことにします（爆）。

実験環境は「Intel Core2 Quad CPU Q8200 2.33GHz」「2GB RAM」で、「VS2008」にて「/O2」でコンパイルしています。テスト画像はかの有名な“lenna”（512×512サイズ、8ビットグレースケール）です。安定した結果で比較できるように5回実行してそれらの中央値を採用しました。結果は以下のとおりでした。なおSNR∞を確認済みです（ただしOpenCVのみ55[dB]前後でした）。

まずは計算時間の比較です。#はコア数です。マルチコア演算の方はおおよそコア数に比例して速くなる感じですね。SSEの方はfloat４並列なんですが２倍程度までしか速くなってません。特にカーネルサイズが小さいと利幅が小さいです（並列化の導入コストが償却できてないのか）。コア数が増えると速度がやや不安定でした（スケジューリングが難くなるからか）。・・・そして小カーネルではOpenCVが倍速orz。きっとオレの実装が悪いんでしょうね（；ε；）。

次にSSEの導入による速度UPの割合を観察します。カーネルが大きくなるに連れて割合がUPしていて、この図の最大ケースで約２倍くらい。もっと大きいカーネルだと理論値である4倍に漸近するのか気になりますが、画像サイズの関係でこれ以上大きいカーネルを試せてませんorz。HDTVサイズくらいのグレースケールのテスト画像、知ってる人プリーズ。

というわけで取り急ぎご報告。早速OpenCVの実装を読まねばならんですな。