C（SIMD）の画像とソーベルフィルタ最適化の高速移転

Question

私は本当に（本当に）高速Sobel operatorを実装したいと思います。私の友人と私は書きました（ソースはhereです）。以下は、私がこれまでに考えていることです...

まず、8ビットの符号なし整数配列にイメージをグレースケールの画像として格納します。

実際のソーベルフィルタを書くには、ピクセルごとにGxとGyを計算する必要があります。これらの数値はそれぞれ、原点の次​​の6ピクセルのおかげで計算されます。しかし、SIMD命令は私に16または32（AVX）ピクセルを扱うことを可能にします。 IによりGyとを計算することができるようにできれば、オペレータのカーネルは、いくつかの素晴らしい特性を有している：それぞれを減算

• IとI + 2行といくつかの他の画像の第i + 1行目（配列で結果を格納その後、追加）
• i + 1及びi + 2列の二回、私を追加することは、私はGxとを計算するために）同じ（ただし、転置を行うだろう

最終画像のI + 1列を与えます2枚の写真。

いくつかの注意：

• すべてが最初に割り当てられますので、私はメモリの割り当てを気にしないでください。
私はオーバーフローに対処し、4によって値を割る問題に署名することができ
• （おかげで_mm_srli_epi8する） (uint8_t >> 2 - uint8_t >> 2) = int7_t //really store as int8_t
int7_t + uint8_t << 1 >> 2 + int7_t = uint8_t
//some precision is lost but I don't care

私が直面してる本当の問題は、行から列へ行くことです。それ以外の場合は、SIMDレジスタに画像をロードできませんでした。私は少なくとも3回イメージを反転する必要がありますか？

元の画像に戻ってください。次に、私はGxとGyの最初のステップを計算してから、結果のピクチャを反転して2番目のステップを計算することができます。

だから、ここに私の質問です：

• は、実装のこの種のは良いアイデアですか？
• ダムアルゴリズムよりも高速に配列を転置する方法はありますか？ （私はそうは思わない）
• ボトルネックはどこですか？ （いずれかの推測？：P）

Answer 1

私は転位/ 2パスはSobelオペレータコードを最適化するのには適していないと思います。 Sobel演算子は計算機能ではないので、転置/ 2パスアクセスのメモリアクセスを無駄にすることはこの場合にはうまくいかない。私はSSEがどのくらい速く得ることができるかを見るためにいくつかのSobel Operatorテストコードを書きました。このコードは最初と最後のエッジピクセルを処理せず、FPUを使用してsqrt（）値を計算します。

ソーベル演算子は、乗算、1平方根、11加算、2分/最大、12リードアクセスおよび1ライトアクセス演算子を必要とする。つまり、コードを最適化すれば、20〜30サイクルでコンポーネントを処理できます。

FloatSobel（）関数は、2113044のCPUサイクルを処理するために256 * 256の画像処理を32.76サイクル/コンポーネントで処理しました。私はこのサンプルコードをSSEに変換します。

void FPUSobel() 
{ 
    BYTE* image_0 = g_image + g_image_width * 0; 
    BYTE* image_1 = g_image + g_image_width * 1; 
    BYTE* image_2 = g_image + g_image_width * 2; 
    DWORD* screen = g_screen + g_screen_width*1; 

    for(int y=1; y<g_image_height-1; ++y) 
    { 
     for(int x=1; x<g_image_width-1; ++x) 
     { 
      float gx = image_0[x-1] * (+1.0f) + 
         image_0[x+1] * (-1.0f) + 
         image_1[x-1] * (+2.0f) + 
         image_1[x+1] * (-2.0f) + 
         image_2[x-1] * (+1.0f) + 
         image_2[x+1] * (-1.0f); 

      float gy = image_0[x-1] * (+1.0f) + 
         image_0[x+0] * (+2.0f) + 
         image_0[x+1] * (+1.0f) + 
         image_2[x-1] * (-1.0f) + 
         image_2[x+0] * (-2.0f) + 
         image_2[x+1] * (-1.0f); 


      int result = (int)min(255.0f, max(0.0f, sqrtf(gx * gx + gy * gy))); 

      screen[x] = 0x01010101 * result; 
     } 
     image_0 += g_image_width; 
     image_1 += g_image_width; 
     image_2 += g_image_width; 
     screen += g_screen_width; 
    } 
} 

SseSobel（）関数は、同じ256 * 256イメージを処理するために613220 CPUサイクルを要しました。 FPUSobel（）よりも9.51サイクル/コンポーネント、3.4倍高速でした。最適化するスペースがいくつかありますが、4ウェイSIMDを使用しているため、4倍より速くはありません。

この関数は、4ピクセルを一度に処理するSoAアプローチを使用しました。 AoSを使用するためにトランスポーズ/シャッフルする必要があるため、SoAはほとんどのアレイまたは画像データでAoSより優れています。そしてSoAは、一般的なCコードをSSEコードに簡単に変更することができます。

void SseSobel() 
{ 
    BYTE* image_0 = g_image + g_image_width * 0; 
    BYTE* image_1 = g_image + g_image_width * 1; 
    BYTE* image_2 = g_image + g_image_width * 2; 
    DWORD* screen = g_screen + g_screen_width*1; 

    __m128 const_p_one = _mm_set1_ps(+1.0f); 
    __m128 const_p_two = _mm_set1_ps(+2.0f); 
    __m128 const_n_one = _mm_set1_ps(-1.0f); 
    __m128 const_n_two = _mm_set1_ps(-2.0f); 

    for(int y=1; y<g_image_height-1; ++y) 
    { 
     for(int x=1; x<g_image_width-1; x+=4) 
     { 
      // load 16 components. (0~6 will be used) 
      __m128i current_0 = _mm_unpacklo_epi8(_mm_loadu_si128((__m128i*)(image_0+x-1)), _mm_setzero_si128()); 
      __m128i current_1 = _mm_unpacklo_epi8(_mm_loadu_si128((__m128i*)(image_1+x-1)), _mm_setzero_si128()); 
      __m128i current_2 = _mm_unpacklo_epi8(_mm_loadu_si128((__m128i*)(image_2+x-1)), _mm_setzero_si128()); 

      // image_00 = { image_0[x-1], image_0[x+0], image_0[x+1], image_0[x+2] } 
      __m128 image_00 = _mm_cvtepi32_ps(_mm_unpacklo_epi16(current_0, _mm_setzero_si128())); 
      // image_01 = { image_0[x+0], image_0[x+1], image_0[x+2], image_0[x+3] } 
      __m128 image_01 = _mm_cvtepi32_ps(_mm_unpacklo_epi16(_mm_srli_si128(current_0, 2), _mm_setzero_si128())); 
      // image_02 = { image_0[x+1], image_0[x+2], image_0[x+3], image_0[x+4] } 
      __m128 image_02 = _mm_cvtepi32_ps(_mm_unpacklo_epi16(_mm_srli_si128(current_0, 4), _mm_setzero_si128())); 
      __m128 image_10 = _mm_cvtepi32_ps(_mm_unpacklo_epi16(current_1, _mm_setzero_si128())); 
      __m128 image_12 = _mm_cvtepi32_ps(_mm_unpacklo_epi16(_mm_srli_si128(current_1, 4), _mm_setzero_si128())); 
      __m128 image_20 = _mm_cvtepi32_ps(_mm_unpacklo_epi16(current_2, _mm_setzero_si128())); 
      __m128 image_21 = _mm_cvtepi32_ps(_mm_unpacklo_epi16(_mm_srli_si128(current_2, 2), _mm_setzero_si128())); 
      __m128 image_22 = _mm_cvtepi32_ps(_mm_unpacklo_epi16(_mm_srli_si128(current_2, 4), _mm_setzero_si128())); 

      __m128 gx = _mm_add_ps(_mm_mul_ps(image_00,const_p_one), 
         _mm_add_ps(_mm_mul_ps(image_02,const_n_one), 
         _mm_add_ps(_mm_mul_ps(image_10,const_p_two), 
         _mm_add_ps(_mm_mul_ps(image_12,const_n_two), 
         _mm_add_ps(_mm_mul_ps(image_20,const_p_one), 
            _mm_mul_ps(image_22,const_n_one)))))); 

      __m128 gy = _mm_add_ps(_mm_mul_ps(image_00,const_p_one), 
         _mm_add_ps(_mm_mul_ps(image_01,const_p_two), 
         _mm_add_ps(_mm_mul_ps(image_02,const_p_one), 
         _mm_add_ps(_mm_mul_ps(image_20,const_n_one), 
         _mm_add_ps(_mm_mul_ps(image_21,const_n_two), 
            _mm_mul_ps(image_22,const_n_one)))))); 

      __m128 result = _mm_min_ps(_mm_set1_ps(255.0f), 
          _mm_max_ps(_mm_set1_ps(0.0f), 
             _mm_sqrt_ps(_mm_add_ps(_mm_mul_ps(gx, gx), _mm_mul_ps(gy,gy))))); 

      __m128i pack_32 = _mm_cvtps_epi32(result); //R32,G32,B32,A32 
      __m128i pack_16 = _mm_packs_epi32(pack_32, pack_32); //R16,G16,B16,A16,R16,G16,B16,A16 
      __m128i pack_8 = _mm_packus_epi16(pack_16, pack_16); //RGBA,RGBA,RGBA,RGBA 
      __m128i unpack_2 = _mm_unpacklo_epi8(pack_8, pack_8); //RRGG,BBAA,RRGG,BBAA 
      __m128i unpack_4 = _mm_unpacklo_epi8(unpack_2, unpack_2); //RRRR,GGGG,BBBB,AAAA 

      _mm_storeu_si128((__m128i*)(screen+x),unpack_4); 
     } 
     image_0 += g_image_width; 
     image_1 += g_image_width; 
     image_2 += g_image_width; 
     screen += g_screen_width; 
    } 
}