NEONを使って画像処理を高速化する (90度回転処理)

ソフトウェア

前回のグレースケール化と一緒に書こうと思っていた90度回転処理の高速化について書きます。

前回の話はこちら:

画像の90度回転処理について考える

今回とりあげるのは画像の90度回転処理です。それをNEONを使用して高速化させるテクニックについて書きます。

画像の回転処理はいろいろあり、任意の角度に回転させる、任意の座標を中心にして回転させる、といったことも考えられますが今回は簡単な場合だけを考えます。

今回対象とする処理は時計回りに90度だけ回転させる処理です。回転させると左側が上になり右側が下になるだけのとても簡単な処理のみを対象とします。下図に例を載せます。

この回転を行う処理をC言語で単純に実装すると次のようになります。

#include <stddef.h>

void rotate90(unsigned char *src, unsigned char *dst, size_t width, size_t height)
{
        for(size_t h = 0; h < height; h++) {
                for(size_t w = 0; w < width; w++) {
                        dst[(w * width + (height - h - 1)) * 3 + 0] = src[(h * width + w) * 3 + 0];
                        dst[(w * width + (height - h - 1)) * 3 + 1] = src[(h * width + w) * 3 + 1];
                        dst[(w * width + (height - h - 1)) * 3 + 2] = src[(h * width + w) * 3 + 2];
                }
        }
}

1ピクセルずつR, G, Bの3要素を回転後の座標にコピーするだけの処理です。回転の前後で幅と高さが入れ替わることに注意します。

NEONを使用して高速化する

並列化しやすいデータ処理を考える

次にNEONを使用した高速化について考えます。とは言っても上記のC言語実装の処理をそのまま並列化すると回転前の横のデータはメモリ空間で連続しているため良いものの、回転後の縦のデータになるとメモリ空間では飛び飛びのデータになってしまうため思うように高速化できません。
そこで回転の処理そのものを工夫します。

まずは問題を簡単にするために4×4の画像について考えます。これを時計回りに90度回転すると次のようになります。ここでアルファベットの1文字が1ピクセルを表すものとします。

A B C D       M I E A
E F G H  ==>  N J F B
I J K L  ==>  O K G C
M N O P       P L H D

この回転と同じ結果になるように次のように処理します。まずは4×4の画像を2×2の画像4つと捉えて、2×2の画像を行列として考え転置をとります。具体的には次のようになります。

まずは(考える上で)4×4の画像を2×2の画像4つに分解します。

A B | C D
E F | G H
- - + - -
I J | K L
M N | O P

次にこの4つの行列に対して転置をとります。
上が転置前の2×2画像4つを表していて下が転置後です。太字の文字は転置で移動したものを表します。

A B | C D | I J | K L
E F | G H | M N | O P
- - + - - + - - + - -
A E | C G | I M | K O
B F | D H | J N | L P

転置をとったら元の4×4の位置に戻します。

A E | C G
B F | D H
- - + - -
I M | K O
J N | L P

次にこの4×4の行列の左右を入れ替えます。同じ行の中で左右(列)のみ入れ替えです。
"A E C G""G C E A"となるようにします。入れ替え後のイメージは次の通りです。

G C E A
H D F B
O K M I
P L N J

最後に1行目の前2列と3行目の後ろ2列を、2行目の前2列と4行目の後ろ2列を入れ替えます。

G C E A | M I E A
H D F B | N J F B
O K M I | O K G C
P L N J | P L H D

ここで入れ替え後の左側の4×4画像を見ると、目的だった90度回転した画像を得ることができました。

NEONの命令で考える

ここまでの回転処理をNEONの命令に当てはめて考えてみます。
まず転置の処理を実現する転置命令がNEONの命令セットにあります。”trn1″命令と”trn2″命令です。転置を意味する”transpose”そのままの名称です。
例としてベクトルレジスタv0とv1に8bitのデータが8個ずつ入っているとします。(4つではないのは8bit x 4 = 32bitだけを扱うベクトルレジスタの操作に対応していないため、64bitならば扱うことができる。)この場合各レジスタから2つずつを取り出していき2×2の行列が4つあるとみなして転置が行われます。具体的には次のようになります。

レジスタv0: 0 1 2 3 4 5 6 7
レジスタv1: 8 9 a b c d e f
- - - - - - - - - - - - - -
trn1 v2.b8, v0.8b, v1.8b
trn2 v3.b8, v0.8b, v1.8b
- - - - - - - - - - - - - -
レジスタv2: 0 8 2 a 4 c 6 e
レジスタv3: 1 9 3 b 5 d 7 f

次に列の前後の入れ替えは”rev64″命令を使用します。具体的には下記のようになります。ここでrev64命令の上が実行前のデータで下が実行後のデータを表します。

レジスタv2: 0 8 2 a 4 c 6 e
- - - - - - - - - - - - - -
rev64 v2.8b, v2.8b
- - - - - - - - - - - - - -
レジスタv2: e 6 c 4 a 2 8 0

最後に2要素ずつの入れ替えは”tbx”命令を使用することで実現できます。tbx命令はルックアップテーブルを参照して要素を取り出ししてくれますので事前に2要素ずつ入れ替えるルックアップテーブルを作成しておきtbx命令を使用します。

レジスタv9: 16 17 18 19  0  1  2  3
レジスタv2:  e  6  c  4  a  2  8  0
レジスタv3:  f  7  d  5  b  3  9  1
- - - - - - - - - - - - - - - - - -
tbx v4.8b, {v2.16b, v3.16b}, v9.8b
v2とv3の結: xx xx xx xx xx xx xx xx  e  6  c  4  a  2  8  0  xx xx xx xx xx xx xx xx  f  7  d  5  b  3  9  1
- - - - - - - - - - - - - - - - - -
レジスタv4:  0  8  2  a  1  9  3  b

ここでレジスタv9が事前に値をセットしておくルックアップテーブルとなります。
レジスタv2とv3は並び替え前のデータでレジスタv4が出力先になります。

注意として、tbx命令はレジスタを2つ結合してそのインデックスを取るためルックアップテーブルの指定に気をつけることと、要素が8つでも8bit x 16個 (16b)しか指定できないため上の”xx”で表した不要なデータまでカウントされることがあげられます。

NEONを利用した実装

それでは実際にアセンブラでNEONの命令を用いてここまでの処理を実装していきます。

	.text
	.align 2
	.p2align 4,,15
	.global rotate90_8x8
rotate90_8x8:
	# x2 = in_width * 3
	# x3 = out_width * 3
	add	x2, x2, x2, lsl 1
	add	x3, x3, x3, lsl 1
	# load look-up table
	adrp	x11, .Ltable1
	ldr	x8, [x11, #:lo12:.Ltable1]
	adrp	x11, .Ltable2
	ldr	x9, [x11, #:lo12:.Ltable2]
	# load above 8x4 pixels
	ld3	{v0.8b - v2.8b}, [x0]
	add	x4, x0, x2
	ld3	{v3.8b - v5.8b}, [x4]
	add	x4, x4, x2
	ld3	{v6.8b - v8.8b}, [x4]
	add	x4, x4, x2
	ld3	{v9.8b - v11.8b}, [x4]
	# load below 8x4 pixels
	add	x4, x4, x2
	ld3	{v12.8b - v14.8b}, [x4]
	add	x4, x4, x2
	ld3	{v15.8b - v17.8b}, [x4]
	add	x4, x4, x2
	ld3	{v18.8b - v20.8b}, [x4]
	add	x4, x4, x2
	ld3	{v21.8b - v23.8b}, [x4]
	# Red (ch1)
	##
	trn1	v24.8b, v0.8b, v3.8b
	trn2	v25.8b, v0.8b, v3.8b
	trn1	v26.8b, v6.8b, v9.8b
	trn2	v27.8b, v6.8b, v9.8b
	trn1	v28.8b, v12.8b, v15.8b
	trn2	v29.8b, v12.8b, v15.8b
	trn1	v30.8b, v18.8b, v21.8b
	trn2	v31.8b, v18.8b, v21.8b
	rev64	v24.8b, v24.8b
	rev64	v25.8b, v25.8b
	rev64	v26.8b, v26.8b
	rev64	v27.8b, v27.8b
	rev64	v28.8b, v28.8b
	rev64	v29.8b, v29.8b
	rev64	v30.8b, v30.8b
	rev64	v31.8b, v31.8b
	ext	v3.8b, v24.8b, v24.8b, #4
	ext	v9.8b, v25.8b, v25.8b, #4
	ext	v0.8b, v26.8b, v26.8b, #4
	ext	v6.8b, v27.8b, v27.8b, #4
	ext	v15.8b, v28.8b, v28.8b, #4
	ext	v21.8b, v29.8b, v29.8b, #4
	ext	v12.8b, v30.8b, v30.8b, #4
	ext	v18.8b, v31.8b, v31.8b, #4
	trn1	v29.4h, v0.4h, v3.4h
	trn2	v25.4h, v0.4h, v3.4h
	trn1	v31.4h, v6.4h, v9.4h
	trn2	v27.4h, v6.4h, v9.4h
	trn1	v28.4h, v12.4h, v15.4h
	trn2	v24.4h, v12.4h, v15.4h
	trn1	v30.4h, v18.4h, v21.4h
	trn2	v26.4h, v18.4h, v21.4h
	# load look-up table
	mov	x6, v24.d[0]
	mov	x7, v25.d[0]
	mov	v24.d[0], x8
	mov	v25.d[0], x9
	# tbx
	tbx	v3.8b, {v26.16b, v27.16b}, v24.8b
	tbx	v15.8b, {v26.16b, v27.16b}, v25.8b
	tbx	v6.8b, {v28.16b, v29.16b}, v24.8b
	tbx	v18.8b, {v28.16b, v29.16b}, v25.8b
	tbx	v9.8b, {v30.16b, v31.16b}, v24.8b
	tbx	v21.8b, {v30.16b, v31.16b}, v25.8b
	mov	v26.d[0], x6
	mov	v27.d[0], x7
	tbx	v0.8b, {v26.16b, v27.16b}, v24.8b
	tbx	v12.8b, {v26.16b, v27.16b}, v25.8b
	# Green (ch2)
	##
	trn1	v24.8b, v1.8b, v4.8b
	trn2	v25.8b, v1.8b, v4.8b
	trn1	v26.8b, v7.8b, v10.8b
	trn2	v27.8b, v7.8b, v10.8b
	trn1	v28.8b, v13.8b, v16.8b
	trn2	v29.8b, v13.8b, v16.8b
	trn1	v30.8b, v19.8b, v22.8b
	trn2	v31.8b, v19.8b, v22.8b
	rev64	v24.8b, v24.8b
	rev64	v25.8b, v25.8b
	rev64	v26.8b, v26.8b
	rev64	v27.8b, v27.8b
	rev64	v28.8b, v28.8b
	rev64	v29.8b, v29.8b
	rev64	v30.8b, v30.8b
	rev64	v31.8b, v31.8b
	ext	v4.8b, v24.8b, v24.8b, #4
	ext	v10.8b, v25.8b, v25.8b, #4
	ext	v1.8b, v26.8b, v26.8b, #4
	ext	v7.8b, v27.8b, v27.8b, #4
	ext	v16.8b, v28.8b, v28.8b, #4
	ext	v22.8b, v29.8b, v29.8b, #4
	ext	v13.8b, v30.8b, v30.8b, #4
	ext	v19.8b, v31.8b, v31.8b, #4
	trn1	v29.4h, v1.4h, v4.4h
	trn2	v25.4h, v1.4h, v4.4h
	trn1	v31.4h, v7.4h, v10.4h
	trn2	v27.4h, v7.4h, v10.4h
	trn1	v28.4h, v13.4h, v16.4h
	trn2	v24.4h, v13.4h, v16.4h
	trn1	v30.4h, v19.4h, v22.4h
	trn2	v26.4h, v19.4h, v22.4h
	# load look-up table
	mov	x6, v24.d[0]
	mov	x7, v25.d[0]
	mov	v24.d[0], x8
	mov	v25.d[0], x9
	# tbx
	tbx	v4.8b, {v26.16b, v27.16b}, v24.8b
	tbx	v16.8b, {v26.16b, v27.16b}, v25.8b
	tbx	v7.8b, {v28.16b, v29.16b}, v24.8b
	tbx	v19.8b, {v28.16b, v29.16b}, v25.8b
	tbx	v10.8b, {v30.16b, v31.16b}, v24.8b
	tbx	v22.8b, {v30.16b, v31.16b}, v25.8b
	mov	v26.d[0], x6
	mov	v27.d[0], x7
	tbx	v1.8b, {v26.16b, v27.16b}, v24.8b
	tbx	v13.8b, {v26.16b, v27.16b}, v25.8b
	# Blue (ch3)
	##
	trn1	v24.8b, v2.8b, v5.8b
	trn2	v25.8b, v2.8b, v5.8b
	trn1	v26.8b, v8.8b, v11.8b
	trn2	v27.8b, v8.8b, v11.8b
	trn1	v28.8b, v14.8b, v17.8b
	trn2	v29.8b, v14.8b, v17.8b
	trn1	v30.8b, v20.8b, v23.8b
	trn2	v31.8b, v20.8b, v23.8b
	rev64	v24.8b, v24.8b
	rev64	v25.8b, v25.8b
	rev64	v26.8b, v26.8b
	rev64	v27.8b, v27.8b
	rev64	v28.8b, v28.8b
	rev64	v29.8b, v29.8b
	rev64	v30.8b, v30.8b
	rev64	v31.8b, v31.8b
	ext	v5.8b, v24.8b, v24.8b, #4
	ext	v11.8b, v25.8b, v25.8b, #4
	ext	v2.8b, v26.8b, v26.8b, #4
	ext	v8.8b, v27.8b, v27.8b, #4
	ext	v17.8b, v28.8b, v28.8b, #4
	ext	v23.8b, v29.8b, v29.8b, #4
	ext	v14.8b, v30.8b, v30.8b, #4
	ext	v20.8b, v31.8b, v31.8b, #4
	trn1	v29.4h, v2.4h, v5.4h
	trn2	v25.4h, v2.4h, v5.4h
	trn1	v31.4h, v8.4h, v11.4h
	trn2	v27.4h, v8.4h, v11.4h
	trn1	v28.4h, v14.4h, v17.4h
	trn2	v24.4h, v14.4h, v17.4h
	trn1	v30.4h, v20.4h, v23.4h
	trn2	v26.4h, v20.4h, v23.4h
	# load look-up table
	mov	x6, v24.d[0]
	mov	x7, v25.d[0]
	mov	v24.d[0], x8
	mov	v25.d[0], x9
	# tbx
	tbx	v5.8b, {v26.16b, v27.16b}, v24.8b
	tbx	v17.8b, {v26.16b, v27.16b}, v25.8b
	tbx	v8.8b, {v28.16b, v29.16b}, v24.8b
	tbx	v20.8b, {v28.16b, v29.16b}, v25.8b
	tbx	v11.8b, {v30.16b, v31.16b}, v24.8b
	tbx	v23.8b, {v30.16b, v31.16b}, v25.8b
	mov	v26.d[0], x6
	mov	v27.d[0], x7
	tbx	v2.8b, {v26.16b, v27.16b}, v24.8b
	tbx	v14.8b, {v26.16b, v27.16b}, v25.8b
	# store
	st3	{v0.8b - v2.8b}, [x1]
	add	x5, x1, x3
	st3	{v3.8b - v5.8b}, [x5]
	add	x5, x5, x3
	st3	{v6.8b - v8.8b}, [x5]
	add	x5, x5, x3
	st3	{v9.8b - v11.8b}, [x5]
	add	x5, x5, x3
	st3	{v12.8b - v14.8b}, [x5]
	add	x5, x5, x3
	st3	{v15.8b - v17.8b}, [x5]
	add	x5, x5, x3
	st3	{v18.8b - v20.8b}, [x5]
	add	x5, x5, x3
	st3	{v21.8b - v23.8b}, [x5]
	ret
.Ldata:
	.section .rodata
	.align	4
.Ltable1:
	.byte 0
	.byte 1
	.byte 2
	.byte 3
	.byte 16
	.byte 17
	.byte 18
	.byte 19
	.align	4
.Ltable2:
	.byte 4
	.byte 5
	.byte 6
	.byte 7
	.byte 20
	.byte 21
	.byte 22
	.byte 23

この実装では変換前データアドレスから8×8ピクセルのRGBデータを読み込んで変換後データアドレスに書き出すようになっています。そのため8×8よりも大きい画像の場合にはこの処理を複数回呼び出す必要があります。また、8ピクセルで割り切れない大きさの場合は余り分を別途処理する必要があります。

処理速度の比較

C言語で実装した回転処理とNEONを利用した回転処理の処理速度を比較してみます。

計測結果は次の通りです。回転に使用した画像のサイズは1920×1080です。それぞれの場合で100回ずつ実行して平均などを集計しました。時間の単位はミリ秒です。

実装方法最適化のレベル平均中央値標準偏差
C言語O254.46554.4250.145
C言語Ofast54.45854.4300.157
アセンブラ43.66843.6200.194
1920×1080の画像を回転

小さい画像として640×480の画像でも試してみました。

実装方法最適化のレベル平均中央値標準偏差
C言語O25.5225.5200.111
C言語Ofast5.5055.4750.141
アセンブラ5.8505.8500.055
640×480の画像を回転

1920×1080のそこそこ大きな画像だとNEONで実装したものが速いという結果になりましたが、640×480と小さめの画像では逆に遅くなりました。

付録

計測時に使用したコードなどを載せておきます。

Makefile

CC=gcc
all:
        $(CC) -c rotate90_8x8.s
        $(CC) -c -O2 rotate90_neon.c
        $(CC) -c -Ofast rotate90_c.c
        $(CC) rotate90_8x8.o rotate90_neon.o rotate90_c.o main.c -o main

アセンブリを呼び出すラッパー

#include <stddef.h>

void rotate90_8x8(unsigned char *src, unsigned char *dst, size_t in_width, size_t out_width);
void rotate90_neon(unsigned char* __restrict src, unsigned char* __restrict dst, size_t width, size_t height)
{
        static const size_t h_step = 8;
        static const size_t w_step = 8;
        for(size_t h = 0; h < height / h_step; h++) {
                for(size_t w = 0; w < width / w_step; w++) {
                        size_t src_index = 3 * ((h * h_step) * width + (w * w_step));
                        size_t dst_index = 3 * ((w * w_step) * height + (height - (h * h_step) - h_step));
                        rotate90_8x8(&src[src_index], &dst[dst_index], width, height);
                }
        }
}

C言語のループ実装 (比較対象)

#include <stddef.h>

void rotate90_c(unsigned char* __restrict src, unsigned char* __restrict dst, size_t width, size_t height)
{
        for(size_t h = 0; h < height; h++) {
                for(size_t w = 0; w < width; w++) {
                        dst[(w * width + (height - h - 1)) * 3 + 0] = src[(h * width + w) * 3 + 0];
                        dst[(w * width + (height - h - 1)) * 3 + 1] = src[(h * width + w) * 3 + 1];
                        dst[(w * width + (height - h - 1)) * 3 + 2] = src[(h * width + w) * 3 + 2];
                }
        }
}

メインルーチン

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

void rotate90_neon(unsigned char* __restrict src, unsigned char* __restrict dst, size_t width, size_t height);
void rotate90_c(unsigned char* __restrict src, unsigned char* __restrict dst, size_t width, size_t height);

int main(int argc, char *argv[])
{
        static const int height = 1920;
        static const int width = 1080;

        unsigned char *src = (unsigned char *)malloc(width * height * 3);
        unsigned char *dst1 = (unsigned char *)malloc(width * height * 3);
        unsigned char *dst2 = (unsigned char *)malloc(width * height * 3);
        struct timespec ts_start, ts_end;

        for(int h = 0; h < height; h++) {
                for(int w = 0; w < width; w++) {
                        src[3 * (h * width + w) + 0] = (255 / width) * w + h;
                        src[3 * (h * width + w) + 1] = (255 / width) * w + h + 1;
                        src[3 * (h * width + w) + 2] = (255 / width) * w + h + 2;
                }
        }

        timespec_get(&ts_start, TIME_UTC);
        rotate90_c(src, dst1, width, height);
        timespec_get(&ts_end, TIME_UTC);
        printf("time: %.2lf\n",
                ((double)(ts_end.tv_sec - ts_start.tv_sec) * 1000 +
                ((double)(ts_end.tv_nsec - ts_start.tv_nsec)) * 0.001 * 0.001));

        timespec_get(&ts_start, TIME_UTC);
        rotate90_neon(src, dst2, width, height);
        timespec_get(&ts_end, TIME_UTC);
        printf("time: %.2lf\n",
                ((double)(ts_end.tv_sec - ts_start.tv_sec) * 1000 +
                ((double)(ts_end.tv_nsec - ts_start.tv_nsec)) * 0.001 * 0.001));

        free(src);
        free(dst1);
        free(dst2);
        return 0;
}

コメント