この記事はOpenCV Advent Calendar 2021の14日目の記事です。
他の記事は目次にまとめられています。

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

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

過去のアドベントカレンダーでも、@tomoaki_teshima 先生が紹介記事を書いてくださっています。 Universal Intrinsic の紹介

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

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

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

■ テスト環境

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

alalekさんのコメントにもある通り、simd_basic.cpp を見るのが最初になります。

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

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

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

#include <opencv2/core/simd_intrinsics.hpp>

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

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

#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

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

サンプルコード（クリックすると開閉します）

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()

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);

}

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

サンプルコードで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のドキュメントが全然整理されていないのが、非常に厳しいです。

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幅を半分に縮小して代入する関数です。

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} を抽出できます。

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の要素を取得する関数です。

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は、指定された番号の要素で、他の要素を全部埋める

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は、構成要素の中の最大値/最長値を返します。

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は、構成要素の合計値を返します。

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 対応状況を簡単に確認してみる。

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された。しかし、バージョン間での動作は非互換...（なんたることか）

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