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

こんにてぃーは。前回の続きです。今回は、大カーネルのときにSSE並列のコードが遅くなる原因について検証し、その改善を図ります。
まずは処理速度低下の原因を明らかにするために処理時間の内訳を観察しましょう。SSE&シングルコアにおける横方向(H)および縦方向(V)の畳み込みと、ついでに画像の転置処理2回分の処理時間を示したのが以下のグラフです。

カーネルのときは(意外にも)縦フィルタのほうが2倍ほど高速な一方、大カーネルでは圧倒的に遅くなっています。並列化してないコードの方では横フィルタ→転置→横フィルタ→転置という手順を踏んでいたために、この縦フィルタの悪影響を受けなかったわけですね(図に示すように転置は当然ながら一定時間で処理できる)。縦と横での結果の違いから推察するに、この速度の急激な低下の原因はメモリアクセスにあると考えられます。CPUキャッシュの仕組み的には、縦フィルタでの飛び飛びのアドレスへのアクセスよりも、横フィルタのような連続したアドレスへのアクセスのほうがベターなはずです。
これを改善するために、ループアンローリングというテクを導入しましょう。これはループ処理を部分的に展開して高速化を図る手法です(実はOpenCV2.3.1の実装でも使われています)。例えば今回のケースの場合、係数との乗算結果をSSEレジスタに累積していきますが、4画素ではなく8画素・16画素ごとに連続した値を処理したほうがキャッシュヒットしやすそうです。幸いワークスペースであるSSEレジスタがまだ余ってるので、そこを利用しましょう。例えば16画素ステップにする場合、縦フィルタのコードは次のようになります。

void convSseV16(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,s0,s1,s2,s3;
        float* p=src+j*w;
        float* q=dst+j*w;
        for(int i=0;i<w;i+=16) //16 step
        {
            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);
            s1=_mm_load_ps(p+i+ 4); s1=_mm_mul_ps(s1,c);
            s2=_mm_load_ps(p+i+ 8); s2=_mm_mul_ps(s2,c);
            s3=_mm_load_ps(p+i+12); s3=_mm_mul_ps(s3,c);
            for(int v=1;v<=r;v++)
            {
                c=_mm_load_ss(&kernel[v]);
                c=_mm_shuffle_ps(c,c,0);
                //横方向に一気に4x4画素処理しちゃう
                x0=_mm_add_ps(_mm_load_ps(p-v*w+i   ),_mm_load_ps(p+v*w+i   )); x0=_mm_mul_ps(x0,c); s0=_mm_add_ps(s0,x0);
                x0=_mm_add_ps(_mm_load_ps(p-v*w+i+ 4),_mm_load_ps(p+v*w+i+ 4)); x0=_mm_mul_ps(x0,c); s1=_mm_add_ps(s1,x0);
                x0=_mm_add_ps(_mm_load_ps(p-v*w+i+ 8),_mm_load_ps(p+v*w+i+ 8)); x0=_mm_mul_ps(x0,c); s2=_mm_add_ps(s2,x0);
                x0=_mm_add_ps(_mm_load_ps(p-v*w+i+12),_mm_load_ps(p+v*w+i+12)); x0=_mm_mul_ps(x0,c); s3=_mm_add_ps(s3,x0);
            }
            _mm_store_ps(q+i   ,s0);
            _mm_store_ps(q+i+ 4,s1);
            _mm_store_ps(q+i+ 8,s2);
            _mm_store_ps(q+i+12,s3);
        }
        for(int i=w-16;i<w;i+=4) //4 step for the rest
        {
            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_add_ps(_mm_load_ps(p-v*w+i),_mm_load_ps(p+v*w+i)); x0=_mm_mul_ps(x0,c); s0=_mm_add_ps(s0,x0);
            }
            _mm_store_ps(q+i,s0);
        }
    }
    delete kernel;
}

ついでに横フィルタもループアンロールしましょう。横フィルタでもループを展開することでガウスカーネルの係数のロード回数を削減できます。

void convSseH16(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);

    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,s1,s2,s3;
        float* p=src+j*w;
        float* q=dst+j*w;
        for(int i=r4;i<=w-r4-16;i+=16) //16画素ステップで処理
        {
            s0=_mm_setzero_ps();
            s1=_mm_setzero_ps();
            s2=_mm_setzero_ps();
            s3=_mm_setzero_ps();
            for(int u=-r;u<=+r;u++)
            {
                //16画素に対し係数のロードが1回で済むので早い
                c=_mm_load_ss(&kernel[r+u]);
                c=_mm_shuffle_ps(c,c,0);
                //横方向に一気に4x4画素処理しちゃう
                x0=_mm_loadu_ps(p+i+u   ); x0=_mm_mul_ps(x0,c); s0=_mm_add_ps(s0,x0);
                x0=_mm_loadu_ps(p+i+u+ 4); x0=_mm_mul_ps(x0,c); s1=_mm_add_ps(s1,x0);
                x0=_mm_loadu_ps(p+i+u+ 8); x0=_mm_mul_ps(x0,c); s2=_mm_add_ps(s2,x0);
                x0=_mm_loadu_ps(p+i+u+12); x0=_mm_mul_ps(x0,c); s3=_mm_add_ps(s3,x0);
            }
            _mm_store_ps(q+i   ,s0);
            _mm_store_ps(q+i+ 4,s1);
            _mm_store_ps(q+i+ 8,s2);
            _mm_store_ps(q+i+12,s3);
        }
        for(int i=w-r4-16;i<w-r4;i+=4) //端数を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;
}

では実験です。実験条件は前回と同じです。以下の図にて、SSE&シングルコアでの4・8・16画素ステップ処理の結果と、最終的なガウシアンフィルタ実装の結果を比較します。


まず上の図ですが、16画素ステップ処理により縦横共に約1.7倍ほど速くなっており、縦フィルタの性能劣化のタイミングも明らかに遅くなっています。ループを展開しただけとは思えない顕著な改善ですね。次に下の図ですが、このSSE部の改良とOpenMP併用の結果、OpenCVに対して2.5倍程度高速となりました。素晴らしい。シングルコアの場合でも全体を通して10%〜20%の性能向上が見られます。相変わらず大カーネル時には、普通の畳込み実装に肩を並べてしまいますが、転置方式にすればもっと良い結果が得られるでしょう(実験は割愛)。ただ個人的には、頻繁に大カーネルを扱うのであれば、根本的にアルゴリズム的な改良を施すほうが吉だとは思います。これについてはまた気が向いたらポストします。
というわけで、ガウシアンフィルタの畳込み演算を並列処理する実験について連続ポストしてみました。並列演算ハードウェアを駆使すれば結構なサイズのカーネルでもゴリ押しでなんとかなるようです。あるパネルディスカッション?でGoogleの人が「アルゴリズムが並列で処理できさえすればなんとかなる」と言ってましたが、それもうなずける結果です。GPUなどと違って、SSEやマルチコアCPUは現代のほぼ全てのパソコンに搭載されているというのも大きな利点ですね。「あと数倍〜10倍弱高速化できたらなぁ・・・」というときに非常に便利な選択肢だと思います。みなさんも是非活用してみてください。