Universal Intrinsic再入門(別角度から)

この記事は[OpenCV Advent Calendar 2021](https://qiita.com/advent-calendar/2021/opencv)の14日目の記事です。
他の記事は目次にまとめられています。


■ Universal Intrinsic再入門(別角度から)
===========================================

OpenCVには、Universal IntrinsicというSIMD実装を容易にするための仕組みがあります。

過去のアドベントカレンダーでも、@tomoaki_teshima 先生が紹介記事を書いてくださっています。 [Universal Intrinsic の紹介](https://qiita.com/tomoaki_teshima/items/8195c25a3f0448ce300b)

Universal Intrinsicを使う事で、コンパイル先のアーキテクチャを気にする事なく[^1]、SIMD実装ができるメリットがあります。
ただ、惜しむべくはUniversal Intrinsicをアプリケーションで使うサンプルがあまり存在していません。

そこで今回は、このUniversal Intrinsic機能をユーザーアプリケーションから使う方法などを紹介していきたいです！
これによって、ちょっとでもユーザーアプリからでも使われるようになるといいなあ・・・と。

[^1]: いや、性能を出すためには、ちゃんと気にしないとダメですけどね！

:::note warn
Universal Intrinsic機能は、OpenCV modulesのための内部的な機能であり、アプリケーションが使うのはtricky、というコメントが以前issue( https://github.com/opencv/opencv/issues/16732 )でありました。
:::


■ テスト環境
=============

![Ubuntu-21.10](https://img.shields.io/badge/Ubuntu-21.10-brightgreen) 　![OpenCV-4.5.4](https://img.shields.io/badge/OpenCV-4.5.4-brightgreen) 


■ どうすればUniversal Intrinsicをアプリケーションから使えるのか？
==================================================================

[alalekさんのコメント](https://github.com/opencv/opencv/issues/16732)にもある通り、[simd_basic.cpp](https://github.com/opencv/opencv/blob/4.x/samples/cpp/simd_basic.cpp) を見るのが最初になります。

OpenCV 4.5.4でコンパイルができる環境であれば、今すぐに使える環境は整っています！

1. `#include <opencv2/core/simd_intrinsics.hpp>`
2. あとは普通に、coreをリンク対象にする

この１行が、今回の記事のメインであり、本体ですね。

```c++
#include <opencv2/core/simd_intrinsics.hpp>
```

◯SIMD非対応環境のことも考える
-------------------------------

すべてのコンパイルするターゲット環境で、SIMDが使えるとは限りません。
マクロ `CV_SIMD` を使ってサポート状況を確認するようにしましょう。
（なお、もっと厳密にはビット幅のマクロもチェックするべきですね・・・）

```c++
#ifdef CV_SIMD
#ifdef CV_SIMD128
    // Factor
    v_uint8x16 a[3];

    a[0] = v_setall_u8( fB );
    a[1] = v_setall_u8( fG );
    a[2] = v_setall_u8( fR );
#endif
#endif
```



下記サンプルコードはパフォーマンス度外視なのでいぢめないでください！！！

<details><summary>サンプルコード（クリックすると開閉します）</summary><div>

```CMakeLists.txt
cmake_minimum_required(VERSION 2.8)
project(test)

find_package(OpenCV REQUIRED)
if(OpenCV_FOUND)
    add_executable( a.out  main.cpp )
    set(CMAKE_CXX_FLAGS "-O2 -std=c++11 -Wall -pg -g")

    target_include_directories(a.out PUBLIC ${OpenCV_INCLUDE_DIRS})
    target_link_libraries(a.out ${OpenCV_LIBS})
endif()
```

```cpp:main.cpp

#define V_TEST

#include <iostream>
#include <opencv2/opencv.hpp>

#ifdef V_TEST
#  include <opencv2/core/simd_intrinsics.hpp>
#endif

using namespace cv;

int main()
{
    cv::Mat img = imread("opencv.jpg");
    uint8_t fB = 18;
    uint8_t fR = 54;
    uint8_t fG = 255 - fB - fR;

#ifdef V_TEST
#ifdef CV_SIMD
#ifdef CV_SIMD128
    // Factor
    v_uint8x16 a[3];

    a[0] = v_setall_u8( fB );
    a[1] = v_setall_u8( fG );
    a[2] = v_setall_u8( fR );
#endif
#endif
#endif

    for ( int y = 0 ; y < img.rows ; y++ )
    {
        int x;
        int offset = ( y * img.cols ) * 3;

        x = 0;
#ifdef V_TEST
#ifdef CV_SIMD
#ifdef CV_SIMD128
        for ( ; x <= img.cols - 16 ; x+=16 )
        {
            unsigned char *ptr = img.data + offset + x * 3;

            // Load memory as uint8x16
            v_uint8x16 b8x16, g8x16, r8x16;
            v_load_deinterleave( ptr, b8x16, g8x16, r8x16 );

            // Expand from uint8x16 to (uint16x8, uint16x8)
            // And multipy factor.
            v_uint16x8 b16x8[2], g16x8[2], r16x8[2];
            v_mul_expand( b8x16, a[0], b16x8[0], b16x8[1] );
            v_mul_expand( g8x16, a[1], g16x8[0], g16x8[1] );
            v_mul_expand( r8x16, a[2], r16x8[0], r16x8[1] );

            // Add sum and shift
            v_uint16x8 t[2];
            t[0] = ( b16x8[0] + g16x8[0] + r16x8[0] ) >> 8;
            t[1] = ( b16x8[1] + g16x8[1] + r16x8[1] ) >> 8;

            // Convert from (uint16x8 , uint16x8) to (uint8x16)
            b8x16 = v_pack( t[0], t[1] );

            // Store memory
            v_store_interleave( ptr, b8x16, b8x16, b8x16 );
        }
#endif
#endif
#endif
        for ( ; x < img.cols; x++ )
        {
            unsigned char *ptr = img.data + offset + x * 3;
            uint16_t t = ( ptr[0] * fB + ptr[1] * fG + ptr[2] * fR ) >> 8;
            ptr[2] = ptr[1] = ptr[0] = (int8_t) t & 0xFF;
        }
    }

    imwrite("opencv2.jpg", img);

}
```

</div></details>

コンパイル結果を見ても、ちゃんとそれなりの命令を使っていますね…

```c:サンプルコードでV_TESTを定義した場合
  402350:       48 8b 9d 48 ff ff ff    mov    -0xb8(%rbp),%rbx
  402357:       f3 0f 6f 7c 3b e0       movdqu -0x20(%rbx,%rdi,1),%xmm7
  40235d:       f3 0f 6f 64 3b f0       movdqu -0x10(%rbx,%rdi,1),%xmm4
  402363:       f3 0f 6f 0c 3b          movdqu (%rbx,%rdi,1),%xmm1
  402368:       66 0f 70 ec ee          pshufd $0xee,%xmm4,%xmm5
  40236d:       66 0f 70 f7 ee          pshufd $0xee,%xmm7,%xmm6
  402372:       66 0f 60 fd             punpcklbw %xmm5,%xmm7
  402376:       66 0f 60 f1             punpcklbw %xmm1,%xmm6
  40237a:       66 0f 70 c9 ee          pshufd $0xee,%xmm1,%xmm1
  40237f:       66 0f 60 e1             punpcklbw %xmm1,%xmm4
  402383:       66 0f 6f ee             movdqa %xmm6,%xmm5
  402387:       66 0f 68 e8             punpckhbw %xmm0,%xmm5

```

■ どんな命令があるのか
=======================

OpenCVでは、このUniversal Intrinsicのドキュメントが全然整理されていないのが、非常に厳しいです。

https://docs.opencv.org/4.x/df/d91/group__core__hal__intrin.html#ga6093cc09443c787193a24ffe4b0b4dcd

`Core functionality » Hardware Acceleration Layer` の中盤あたりに一応コメントがあるので、これに従う。 `Detailed Description` で検索してみてください。

下記は、integer系の、一部を抜粋し、注釈を加えた表になります。

|Operations Types|	uint 8|	int 8|	uint 16|int 16|	uint 32|int 32|comment|
|:---------------|--------|------|---------|------|--------|------|-------|
|load, store|	x|	x|	x|	x|	x|	x|align有り無し|
|　interleave|	x|	x|	x|	x|	x|	x| |
|expand|	x|	x|	x|	x|	x|	x|bit幅2倍拡張 |
|　expand_low|	x|	x|	x|	x|	x|	x|　下位 |
|　expand_high|	x|	x|	x|	x|	x|	x|　上位 |
|　expand_q|	x|	x|	|	|	|	|　4倍拡張|
|add, sub|	x|	x|	x|	x|	x|	x| |
|　add_wrap, sub_wrap|	x|	x|	x|	x|	|	|without saturation|
|mul|	x|	x|	x|	x|	x|	x| |
|　mul_wrap|	x|	x|	x|	x|	|	|without saturation|
|　mul_expand|	x|	x|	x|	x|	x|	|乗算+bit幅2倍拡張|
|compare|	x|	x|	x|	x|	x|	x|比較 |
|shift|	|	|	x|	x|	x|	x| ">>", "<<" |
|dotprod|	|	|	|	x|	|	x| 内積 |
|　dotprod_fast|	|	|	|	x|	|	x| 内積+ fast |
|　dotprod_expand|	x|	x|	x|	x|	|	x| 内積+bit幅2倍拡張 |
|　dotprod_expand_fast|	x|	x|	x|	x|	|	x| 内積+bit幅2倍拡張+fast |
|logical|	x|	x|	x|	x|	x|	x| 対数 |
|min, max|	x|	x|	x|	x|	x|	x| 最大値・最小値 |
|absdiff|	x|	x|	x|	x|	x|	x| 差 |
|　absdiffs|	|	x|	|	x|	|	| 差?|
|reduce|	x|	x|	x|	x|	x|	x| sad/sum/max/min/... |
|mask|	x|	x|	x|	x|	x|	x| |
|pack|	x|	x|	x|	x|	x|	x|bit幅縮小|
|pack_u|	x|	|	x|	|	|	|bit幅縮小|
|　pack_b|	x|	|	|	|	|	| 
|unpack|	x|	x|	x|	x|	x|	x|？？？|
|extract|	x|	x|	x|	x|	x|	x| 要素ずらし|
|extract_n|	x|	x|	x|	x|	x|	x| laneから抽出|
|reverse|	x|	x|	x|	x|	x|	x| 要素の順番を逆順に入れ替え |
|rotate (lanes)|	x|	x|	x|	x|	x|	x| 要素を１つずつrotate |
|cvt_flt32|	|	|	|	|	|	x| |
|cvt_flt64|	|	|	|	|	|	x| |
|transpose4x4|	|	|	|	|	x|	x| |
|broadcast_element|	|	|	|	|	x|	x| 特定要素を全要素にコピー|

この中でちょっとすぐに効果が思いつかなそうなものをいくつか紹介する。

◯expand, pack/pack_u
---------
- expandは、bit幅を2倍に拡張して、それぞれ2つの変数に代入する関数です。
- packは、2要素を引数にして、bit幅を半分に縮小して代入する関数です。

```c++
v_uint8x16 a; // = { A1, A2, A3, A4, A5, A6, A7, A8 } as uint8
v_uint16x8 b1,b2;

v_expand(a,b1,b2);
  // b1 = {A1, A2, A3, A4} as uint16
  // b2 = {A5, A6, A7, A8} as uint16

uint8x16 c;
c = pack(b2,b1); // {A5, A6, A7, A8, A1, A2, A3, A4} as uint8
```

◯extract
--------------
extractは、2要素を引数に、最大lane数だけずらした要素を取得する関数です。
例えば、`{B1,G1,R1,B2,　G2,R2,B3,G3}`, `{R3,G4,B4,R4,　B5,G5,R5,B6}` から、`{B3,G3,R3, B4,G4,R4, B5,G5}` を抽出できます。

```c++
v_int32x4 a; // ={A1,A2,A3,A4}
v_int32x4 b; // ={B1,B2,B3,B4}

v_int32x4 c;
c = v_extract<0>(a,b); // = {A1,A2,A3,A4}
c = v_extract<1>(a,b); // = {A2,A3,A4,B1}
c = v_extract<2>(a,b); // = {A3,A4,B1,B2}
c = v_extract<3>(a,b); // = {A4,B1,B2,B3}
```

◯extract_n
--------------
extractは、1要素を引数に、特定laneの要素を取得する関数です。

```c++
v_int32x4 a; // ={A1,A2,A3,A4}
int r;
r = v_extract_n<0>; // A1
r = v_extract_n<1>; // A2
r = v_extract_n<2>; // A3
r = v_extract_n<3>; // A4
```

◯broadcast_element
--------------
broadcast_elementは、指定された番号の要素で、他の要素を全部埋める

```c++
v_int32x4 a; // ={A1,A2,A3,A4}
v_int32x4 b = v_broardcase_element<0>(a); // = {A1,A1,A1,A1}
v_int32x4 c = v_broardcase_element<1>(a); // = {A2,A2,A2,A2}
v_int32x4 d = v_broardcase_element<2>(a); // = {A3,A3,A3,A3}
v_int32x4 e = v_broardcase_element<3>(a); // = {A4,A4,A4,A4}
```

◯reduce_max/min
--------------
reduce_max/minは、構成要素の中の最大値/最長値を返します。

```c++
v_int32x4 a; // ={A1,A2,A3,A4}
int vmax = v_reduce_max(a); // = max(A1,A2,A3,A4)
int vmin = v_reduce_min(a); // = min(A1,A2,A3,A4)
```

◯reduce_sum
--------------
reduce_sumは、構成要素の合計値を返します。

```c++
v_int32x4 a; // ={A1,A2,A3,A4}
int vsum = v_reduce_sum(a); // = A1+A2+A3+A4
```

◯除算が無いかも…
------------------
加算・減算・乗算はありますが、除算は無いっぽく見えますね。shift命令で置き換えるか、逆数の乗算にする、ですかね。

■ まとめ
=========
OpenCV Universal Intrinsicを、moduleではなくユーザーアプリケーションから使う方法をまとめました。

これで、自分が書いたカーネル関数などもsimd化できますね！！

（今回の私の例みたいに、コンパイラが最適化した方が速いかもしれませんが…）

以上になります。



■ （おまけ）OpenCV Universal IntrinsicのRISC-V対応状況
==================================================================

◯各Extensionでの対応状況。
---------------------------------------------
OpenCVのRISC-V 対応状況を簡単に確認してみる。

https://github.com/opencv/opencv/tree/4.x/modules/core/include/opencv2/core/hal

| Extension | Version | implementation |
|:-|:-|:-|:-|
| V vector Extension | 0.7.1 (Draft)   | intrin_rvv071.hpp    |
| V vector Extension | 1.0             | intrin_rvv.hpp |
| P packed extension | 0.9.10 (Draft)  | N/S               |

なぜ、V vector Extensionが0.7.1と1.0で分かれているのかというと…

- 世間一般で先にDraft版の0.7.1が広まった。
- そのあとで1.0が正式版でReleaseされた。しかし、バージョン間での動作は非互換...（なんたることか）