Edited at
OpenCVDay 5

Universal Intrinsic の紹介


はじめに

これは、OpenCV Advent Calendar 2016 5日目の記事です。関連記事は目次にまとめられています。


Universal Intrinsicの前に



  • SSENEONと呼ばれる命令群(SIMD命令とも呼ぶ)は画像処理にとって非常に重要かつ有効なものであります。

  • 1命令で、同じ演算を複数のピクセルに適用できるため、画像処理との親和性が高いからです。

  • 一方で、SIMD命令はアーキテクチャ、もしくは命令セット固有のものであり、一度書いてしまうと、ポータビリティの弊害となります


  • SIMDを使う上で、演算速度 対 ポータビリティのトレードオフを考慮する必要があります

具体例を見てみましょう


演算速度無視、ポータビリティ重視のコード


  • 以下のコードは、2つの画像を足し合わせてdstに保存するコードです。


normal.cpp

for(int y = 0;y < height;y++)

{
for(int x = 0;x < width;x++)
{
dst.data[y*width+x] = imgA.data[y*width+x] + imgB.data[y*width+x];
}
}

normal_add.png


  • このコードは恐らくどんなアーキテクチャでも動くでしょう。


  • x86でもIA64でもARMでもMIPSでも、アーキテクチャの差異はコンパイラがハンドリングしてくれるので、どのアーキテクチャでも上記コードは動きます。


演算速度重視、ポータビリティ無視のコード


  • さて、ではここでx86SIMD命令である、SSE2を使って、一部高速化してみましょう。


sse.cpp

for(int y = 0;y < height;y++)

{
int offset = y * width;
for(int x = 0;x <= width-8;x+=8)
{
__m128i va = _mm_loadu_si128((const __m128i*)(imgA.data+offset+x)); // ロード1
__m128i vb = _mm_loadu_si128((const __m128i*)(imgB.data+offset+x)); // ロード2
__m128i vd = _mm_adds_epi16(va, vb); // 加算(演算)
_mm_storeu_si128((__m128i*)(dst+offset+x), vd); // ストア
}
}


  • いきなり見づらくなりましたね。。。

  • あまり深入りはしませんが、ループ内の


    • 最初の2行がメモリからレジスタへのロード

    • 3行目が足し算

    • 4行目がレジスタからメモリへのストア



  • となっています。また、画素はshort(16bit幅)としてコードを書いているので、8要素ごとに計算されます。

  • 以下の図では、青と赤がそれぞれimgAimgB、緑がdstを表しています

    simd_add.png


  • 本当は画像端での処理とか、メモリの番地のアライメントを考慮に入れるとか、まだまだ高速化の余地があるのですが、本題からどんどん離れていくので、ここでは割愛します。


  • ここで重要なのは8要素ごとに処理されるので、画像が十分大きければ処理速度(スループット)が8倍高速化される、という点です。


  • 一方で、(こっからが重要)このコードはx86系のコンパイラでしかコンパイルできません。ポータビリティが失われている訳です。


  • 例えば、上司が「チミ、このアプリ、スマホに移植して」とか言い出すとさぁ大変。最近のスマホはARMばかりですからコンパイラを変えるとすぐビルドすらできなくなります。



演算速度重視、ポータビリティも重視のコード


  • さて、そこまでは当然想定される状況な訳で、「じゃあARMでもx86でもどんと来いや!」というコードを予め書くようになる訳です。(ならない?なりましょう)


chaos.cpp

for(int y = 0;y < height;y++)

{
int offset = y * width;
for(int x = 0;x <= width-8;x+=8)
{
#if defined(__x86_64__) || defined(_M_X64) || defined(_M_IX86) || defined(i386)
__m128i va = _mm_loadu_si128((const __m128i*)(imgA.data+offset+x));
__m128i vb = _mm_loadu_si128((const __m128i*)(imgB.data+offset+x));
__m128i vd = _mm_adds_epi16(va, vb);
_mm_storeu_si128((__m128i*)(dst+offset+x), vd);
#elif defined(__arm__) || defined(_M_ARM) || defined(__aarch64__)
int16x8_t va = vld1q_s16(imgA.data+offset+x);
int16x8_t vb = vld1q_s16(imgB.data+offset+x);
int16x8_t vd = vaddq_s16(va, vb);
vst1q(dst.data+offset+x, vd);
#else
dst.data[x+offset] = imgA.data[x+offset] + imgB.data[x+offset];
dst.data[x+offset+1] = imgA.data[x+offset+1] + imgB.data[x+offset+1];
dst.data[x+offset+2] = imgA.data[x+offset+2] + imgB.data[x+offset+2];
dst.data[x+offset+3] = imgA.data[x+offset+3] + imgB.data[x+offset+3];
dst.data[x+offset+4] = imgA.data[x+offset+4] + imgB.data[x+offset+4];
dst.data[x+offset+5] = imgA.data[x+offset+5] + imgB.data[x+offset+5];
dst.data[x+offset+6] = imgA.data[x+offset+6] + imgB.data[x+offset+6];
dst.data[x+offset+7] = imgA.data[x+offset+7] + imgB.data[x+offset+7];
#endif
}
}


  • さらに可読性が落ちた気がしますが、気にしません

  • これでx86でもARMでもビルドできるし、万が一それ以外のアーキテクチャの場合は最後のelseディレクティブのコードがコンパイルされるので、どのアーキテクチャでもサポートされます。

  • 速度もポータビリティも十分なコードですが、そのために可読性が著しく失われた結果となります。


OpenCVの場合


  • さて、OpenCVアドベントカレンダーなので、OpenCVについて説明します。

  • OpenCVは、様々なアーキテクチャで実行されることを想定しています。

  • またいろんな人が読んで、書いて、実行するので、可読性を無視する訳にも行きません。

  • OpenCV 3.0頃までは、SIMD命令はSSEしかまともにサポートされていませんでしたが、OpenCV 3.1 や本稿執筆時のgit版などでは、NEONサポートも積極的に取り入れられています。

  • という訳で、前述の可読性パフォーマンスの両方を担保するような仕組みが必要になります。


  • そこで、OpenCVではUniversal Intrinsicが実装されています。一体どういうものなのか?試しに先程のコードをUniversal Intrinsicで書き換えてみましょう。



universal_intrinsic_example.cpp

for(int y = 0;y < height;y++)

{
int offset = y * width;
for(int x = 0;x <= width-8;x+=8)
{
v_int16x8 va = v_load(imgA.data+offset+x);
v_int16x8 vb = v_load(imgB.data+offset+x);
v_int16x8 vd = va + vb;
v_store(dst.data+offset+x, vd);
}
}


  • どうでしょうか?可読性は格段に上がってると思います(前のがひどすぎた、とも言う)


  • ifdefでの区切りが一切なくなり、またベクタの型がv_int16x8と書いてあると、「short(16bit)が8個のベクトルだな」と分かりやすいし、va + vbという書き方で足し算が行われるのも、非常に分かりやすいと思います。


Universal Intrinsicの中身


  • さて、この処理が結局のところ、どうやって処理されるのか?という点ですが、modules/core/include/opencv2/core/hal を覗いてみると、


    • intrin_neon.hpp

    • intrin_sse.hpp

    • intrin_cpp.hpp



  • の3つのファイルが見えます。


    • ちなみにこのパスに含まれるhalというキーワードは、3.0で華々しくデビューしたにもかかわらず、3.1で要らない子認定された挙句バラされて各モジュール内に組み込まれた、悲しい運命を辿ったHALモジュール(Hardware Acceleration Layer)の名残です。



  • それぞれのファイルの中で、大量のinline関数が宣言されています。例えば、足し算を探してみましょう。まずはintrin_sse.hppから。


intrin_sse.hpp

#define OPENCV_HAL_IMPL_SSE_BIN_OP(bin_op, _Tpvec, intrin) \

inline _Tpvec operator bin_op (const _Tpvec& a, const _Tpvec& b) \
{ \
return _Tpvec(intrin(a.val, b.val)); \
} \
inline _Tpvec& operator bin_op##= (_Tpvec& a, const _Tpvec& b) \
{ \
a.val = intrin(a.val, b.val); \
return a; \
}
// 中略
OPENCV_HAL_IMPL_SSE_BIN_OP(+, v_int16x8, _mm_adds_epi16)


  • こんな部分が見つかります。

  • 非常に見づらいと思いますが、実はマクロでごちゃごちゃしてるので、プリプロセスを通ると、以下のようなコードになります。


universal_intrinsic_example.cpp

inline v_int16x8 operator + (const v_int16x8& a, const v_int16x8& b)

{
return v_int16x8(_mm_adds_epi16(a.val, b.val));
}


  • 単純に2つのベクトルを足し算するinline関数に置き換えられます。(足し算命令は_mm_adds_epi16)

  • さて、一方でARM向けのintrin_neon.hppを見てみましょう


intrin_neon.hpp

#define OPENCV_HAL_IMPL_NEON_BIN_OP(bin_op, _Tpvec, intrin) \

inline _Tpvec operator bin_op (const _Tpvec& a, const _Tpvec& b) \
{ \
return _Tpvec(intrin(a.val, b.val)); \
} \
inline _Tpvec& operator bin_op##= (_Tpvec& a, const _Tpvec& b) \
{ \
a.val = intrin(a.val, b.val); \
return a; \
}
// 中略
OPENCV_HAL_IMPL_NEON_BIN_OP(+, v_int16x8, vqaddq_s16)


  • こいつらも、プリプロセスすると、以下の様に変わります。


universal_intrinsic_example.cpp

inline v_int16x8 operator + (const int16x8& a, const int16x8& b)

{
return v_int16x8(vqaddq_s16(a.val, b.val));
}


  • やはり、16bit幅の足し算、vqaddq_s16が呼ばれています

  • 要は、アーキテクチャの違いをこの時点で吸収している訳です。そのため、


    • 実際のコードではSIMD命令が走るので、速度のパフォーマンスを享受できる

    • ソースコードの可読性はifdefなどを使わずに、可能な限り良い状態で保持できる



  • この両方を実現した訳です。

  • 厳密に言えば、ifdefintrin.hppに押し込めて、コンパイラ/プラットフォームごとにintrin_sse.hppintrin_neon.hppを適切にincludeしている訳です

  • なお、Visual Studioでは基本的にinline展開はReleaseモードでしか行われないので、Debug版では速度に関して享受できない点に注意しましょう。(多分gccでも同じ)


Universal Intrinsicの未来


  • Universal Intrinsic のおかげで可読性と速度は両立できてる訳ですが、一方で新しいアーキテクチャへの対応が容易になる、というメリットも忘れてはいけません。

  • 筆者が把握してる限りでは、MIPSアーキテクチャがSIMD命令をRelease 5 で発表しました12


    • このSIMD命令を使う場合、OpenCV全体に渡って#ifdef MSA みたいなコードを埋め込む必要はなく、intrin_msa.hppみたいなファイルを実装して追加すれば良いわけです。



  • 現に、OpenCVのissue では、PowerPC版のSIMD最適化求む! みたいなのも募集がかかっていたりします。


    • 金銭的なサポートも検討するぜ、って書いてあるので暇があって金がない人は真面目に考えてみても良いと思います




その他テクい話



  • だんでらいおん先生に煽られる前にもうちょっと読んで役に立つ話を書いておこうと思います


  • SIMD命令のSSENEONは同じ128bit幅の命令ながら、演算内容に違いが多々あります

  • その一つが、shuffle、もしくはpermutation と呼ばれるベクトル内の要素の並べ替えの対応です


  • SSENEONに比べるとシャッフル系の命令が多いです

  • このシャッフルが多用されるのが、カラー画像に対するフィルタ処理です


SSEの場合


  • フィルタリングなどは、RGB各チャンネルそれぞれに独立して適用する必要がありますが、各色はインタリーブしてメモリ上に保存されていますので、とあるRの画素の隣のR画素は、3つ隣に存在します

  • 一番簡単なのはsplit関数を使って各成分ごとに、3枚別の画像に切り出すことなのですが、速度を重視してるときに画像まるまるコピーするような遅い処理をしている暇はありません

  • そこで、unpack命令を使って、デインタリーブ処理をします

  • しかし、このunpack命令は基本的に偶数/奇数をデインタリーブするしか無く、RGBの3色をデインタリーブする命令は存在しません

  • 結果として、SSEを使ってRGBの各チャンネルをロードするには、unpack命令を4段階、計24回実行して各チャンネルごとのベクトルとしてロードします
    deinterleave_sse.png

  • 全部で5段の画像にしましたので、4段の処理が必要。各段の処理は6回unpack命令を実行しています


NEONの場合



  • 一方NEONにはシャッフル命令がありません。嘘。シャッフル系の命令は存在する3。しかし、64bit境界をまたいでの操作は難しく、SSEほど便利ではない。

  • では前述のような各要素のデインタリーブをしたい場合はどうするんでしょうか?

  • 実はNEONのロード命令にはVLD1からVLD4まで4種類存在し、それを使うことで一発でデインタリーブしながらレジスタにロードできます4


  • VLD1はデインタリーブ無しですが、VLD2は要素を2つずつ、VLD3は要素を3つずつ、VLD4は要素を4つずつ、それぞれデインタリーブしながらロードしてくれます


  • VLD3を設計した人とか、神なんじゃないかと思えますね
    deinterleave_neon.png

  • このあたりがコード量の差で見えてきます


intrin_sse.hpp

inline void v_load_deinterleave(const uchar* ptr, v_uint8x16& a, v_uint8x16& b, v_uint8x16& c)

{
__m128i t00 = _mm_loadu_si128((const __m128i*)ptr);
__m128i t01 = _mm_loadu_si128((const __m128i*)(ptr + 16));
__m128i t02 = _mm_loadu_si128((const __m128i*)(ptr + 32));

__m128i t10 = _mm_unpacklo_epi8(t00, _mm_unpackhi_epi64(t01, t01));
__m128i t11 = _mm_unpacklo_epi8(_mm_unpackhi_epi64(t00, t00), t02);
__m128i t12 = _mm_unpacklo_epi8(t01, _mm_unpackhi_epi64(t02, t02));

__m128i t20 = _mm_unpacklo_epi8(t10, _mm_unpackhi_epi64(t11, t11));
__m128i t21 = _mm_unpacklo_epi8(_mm_unpackhi_epi64(t10, t10), t12);
__m128i t22 = _mm_unpacklo_epi8(t11, _mm_unpackhi_epi64(t12, t12));

__m128i t30 = _mm_unpacklo_epi8(t20, _mm_unpackhi_epi64(t21, t21));
__m128i t31 = _mm_unpacklo_epi8(_mm_unpackhi_epi64(t20, t20), t22);
__m128i t32 = _mm_unpacklo_epi8(t21, _mm_unpackhi_epi64(t22, t22));

a.val = _mm_unpacklo_epi8(t30, _mm_unpackhi_epi64(t31, t31));
b.val = _mm_unpacklo_epi8(_mm_unpackhi_epi64(t30, t30), t32);
c.val = _mm_unpacklo_epi8(t31, _mm_unpackhi_epi64(t32, t32));
}



  • 一方ロシアは鉛筆を使ったNEON


intrin_neon.hpp

inline void v_load_deinterleave(const uchar* ptr, v_uint8x16& a, v_uint8x16& b, v_uint8x16& c)

{
uint8x16x3_t v = vld3q_u8(ptr);
a.val = v.val[0];
b.val = v.val[1];
c.val = v.val[2];
}


  • これは便利!


  • NEON版の最後の3行はOpenCVのUniversal Intrinsic構造体に書き戻すための処理ですので、実際の処理はSSE版が15行なのに対し、NEON版では1行で済んでいます


まとめ


  • Universal Intrinsic を紹介して、可読性を保ちながらSIMDで高速化する方法を紹介しました。

  • OpenCVは広く様々なアーキテクチャで使われることを想定したライブラリなので、Universal Intrinsic という妥協点に落ち着いたのでしょう

  • 明日も私の担当で、記事の内容は「ハミング距離の計算はホントに速いのか?(ARM版)」です


備考

筆者は以下の環境で試しました


  • OpenCV masterバージョン (8213e57f)

  • Visual Studio 2012 Update 5

  • Windows 7 Ultimate 64bit