Universal Intrinsic の紹介

# はじめに
これは、OpenCV Advent Calendar 2016 5日目の記事です。関連記事は[目次](http://qiita.com/advent-calendar/2016/opencv)にまとめられています。

# Universal Intrinsicの前に
- ```SSE```や```NEON```と呼ばれる命令群(```SIMD```命令とも呼ぶ)は画像処理にとって非常に重要かつ有効なものであります。
- 1命令で、同じ演算を複数のピクセルに適用できるため、画像処理との親和性が高いからです。
- 一方で、```SIMD```命令はアーキテクチャ、もしくは命令セット固有のものであり、一度書いてしまうと、ポータビリティの弊害となります
- ```SIMD```を使う上で、**演算速度 対 ポータビリティのトレードオフ**を考慮する必要があります

具体例を見てみましょう

## 演算速度無視、ポータビリティ重視のコード
- 以下のコードは、2つの画像を足し合わせて```dst```に保存するコードです。

```cpp: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](https://qiita-image-store.s3.amazonaws.com/0/12162/f1a092ec-41ee-1134-f078-4998b53f1ae8.png "通常の加算処理")

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

## 演算速度重視、ポータビリティ無視のコード
- さて、ではここで```x86```の```SIMD```命令である、```SSE2```を使って、一部高速化してみましょう。

```cpp: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要素ごとに計算されます。
- 以下の図では、青と赤がそれぞれ```imgA```と```imgB```、緑が```dst```を表しています
![simd_add.png](https://qiita-image-store.s3.amazonaws.com/0/12162/c6cf4cb5-5c15-6c46-dfff-bff941ae04d2.png "SIMD命令で加算")

- 本当は画像端での処理とか、メモリの番地のアライメントを考慮に入れるとか、まだまだ高速化の余地があるのですが、本題からどんどん離れていくので、ここでは割愛します。
- ここで重要なのは**8要素ごとに処理される**ので、画像が十分大きければ処理速度(スループット)が**8倍高速化される**、という点です。
- 一方で、（こっからが重要）このコードは```x86```系のコンパイラでしかコンパイルできません。ポータビリティが失われている訳です。
- 例えば、上司が「**チミ、このアプリ、スマホに移植して**」とか言い出すとさぁ大変。最近のスマホは```ARM```ばかりですからコンパイラを変えるとすぐビルドすらできなくなります。

## 演算速度重視、ポータビリティも重視のコード
- さて、そこまでは当然想定される状況な訳で、「じゃあ```ARM```でも```x86```でもどんと来いや！」というコードを予め書くようになる訳です。（ならない？なりましょう）

```cpp: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で書き換えてみましょう。

```cpp: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```から。

```cpp: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)
```

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

```cpp: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```を見てみましょう

```cpp: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)
```

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

```cpp: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```などを使わずに、可能な限り良い状態で保持できる
- この両方を実現した訳です。
- 厳密に言えば、```ifdef``` を```intrin.hpp```に押し込めて、コンパイラ/プラットフォームごとに```intrin_sse.hpp```か```intrin_neon.hpp```を適切に```include```している訳です
- なお、Visual Studioでは基本的にinline展開はReleaseモードでしか行われないので、Debug版では速度に関して享受できない点に注意しましょう。（多分gccでも同じ）

## Universal Intrinsicの未来
- Universal Intrinsic のおかげで可読性と速度は両立できてる訳ですが、一方で新しいアーキテクチャへの対応が容易になる、というメリットも忘れてはいけません。
- 筆者が把握してる限りでは、MIPSアーキテクチャがSIMD命令をRelease 5 で発表しました[^1][^2]。
 - このSIMD命令を使う場合、OpenCV全体に渡って```#ifdef MSA``` みたいなコードを埋め込む必要はなく、```intrin_msa.hpp```みたいなファイルを実装して追加すれば良いわけです。
- 現に、OpenCVのissue では、[PowerPC版のSIMD最適化求む！](https://github.com/opencv/opencv/issues/7207) みたいなのも募集がかかっていたりします。 
 - 金銭的なサポートも検討するぜ、って書いてあるので暇があって金がない人は真面目に考えてみても良いと思います
 
[^1]: [MIPS SIMD architecture - Latest Release of MIPS Architecture Includes Virtualization and SIMD](https://imgtec.com/news/press-release/imagination-reveals-first-mips-warrior-p-class-cpu-core/)
[^2]: [Imagination shares positive MIPS progress - Imagination Technologies] (https://imgtec.com/news/press-release/imagination-shares-positive-mips-progress/)


# その他テクい話
- ~~[だんでらいおん先生]( https://twitter.com/dandelion1124/ )に煽られる前に~~もうちょっと読んで役に立つ話を書いておこうと思います
- ```SIMD```命令の```SSE```と```NEON```は同じ128bit幅の命令ながら、演算内容に違いが多々あります
- その一つが、shuffle、もしくはpermutation と呼ばれる**ベクトル内の要素の並べ替え**の対応です
- ```SSE```は```NEON```に比べるとシャッフル系の命令が多いです
- このシャッフルが多用されるのが、**カラー画像に対するフィルタ処理**です

## SSEの場合
- フィルタリングなどは、RGB各チャンネルそれぞれに独立して適用する必要がありますが、各色はインタリーブしてメモリ上に保存されていますので、とある**Rの画素の隣のR画素**は、3つ隣に存在します
- 一番簡単なのは```split```関数を使って各成分ごとに、3枚別の画像に切り出すことなのですが、速度を重視してるときに画像まるまるコピーするような**遅い**処理をしている暇はありません
- そこで、```unpack```命令を使って、デインタリーブ処理をします
- しかし、この```unpack```命令は基本的に偶数/奇数をデインタリーブするしか無く、RGBの**3色**をデインタリーブする命令は存在しません
- 結果として、```SSE```を使ってRGBの各チャンネルをロードするには、```unpack```命令を4段階、計24回実行して各チャンネルごとのベクトルとしてロードします
![deinterleave_sse.png](https://qiita-image-store.s3.amazonaws.com/0/12162/1cbcd8a9-dc77-dab2-94cb-4a44b286f700.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](https://qiita-image-store.s3.amazonaws.com/0/12162/eedb874a-582c-3650-6598-c2f9343d1cf0.png "神！")
- このあたりがコード量の差で見えてきます

```cpp: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));
}
```
[^3]: [ARM Processors: Coding for NEON - Part 5: Rearr... | ARM Connected Community](https://community.arm.com/groups/processors/blog/2012/03/13/coding-for-neon--part-5-rearranging-vectors)
[^4]: [ARM Processors: Coding for NEON - Part 1: Load ... | ARM Connected Community](https://community.arm.com/groups/processors/blog/2010/03/17/coding-for-neon--part-1-load-and-stores)
- 一方~~ロシアは鉛筆を使った~~```NEON```は

```cpp: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版)](http://qiita.com/tomoaki_teshima/items/00ab8036195b445d1dde)」です

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

- OpenCV masterバージョン (8213e57f)
- Visual Studio 2012 Update 5
- Windows 7 Ultimate 64bit