画像処理を高速化させたい
以前ARM系プロセッサーのSIMDであるNEONを使用して処理を並列化、高速化させることを行いました。今回はARMプロセッサー上で画像処理をする機会があり単純な実装だと遅いと感じたのでNEONを使用した高速化を試みました。
色空間の変換
まず、今回行いたい処理は色空間の変換(RGBから輝度Yへ変換)と時計回りに90度の回転です。この2つの処理に絞ってチューニングしてより処理速度を速くさせるようにします。
初めに色空間の変換の処理を考えます。変換自体は単純で1ピクセルごとにRGBの値に係数をかけ輝度Yに変換します。具体的には下記の式のように変換します。
Y = 0.299 * R + 0.587 * G + 0.114 * BここでYは輝度、R、G、Bがそれぞれ赤緑青の値です。また、全ての値域は0から254です。
通常の実装
この変換をC言語で単純に実装すると次のようになるかと思います。
void convertColorSpaceRGBToY(unsigned char *src, unsigned char *dst, unsigned int len)
{
for(unsigned int i = 0; i < len; i++) {
dst[i] = src[i * 3] * 0.299f + src[i * 3 + 1] * 0.587f + src[i * 3 + 2] * 0.114f;
}
}動きとしてはピクセル数lenを受け取ってその分だけ先程の変換式の通りに変換していきます。変換前、変換後ともにunsigned char型ですが変換途中の計算ではfloat型にキャストして計算しています。
この処理の流れを考えるとまずRGBの値をメモリから3バイト読み取りそれぞれをfloat型に変換、係数をかけて合計を求め、今後はunsigned charに変換し1バイトの値をメモリに書き込みます。その後ループの継続判定としてカウンタをインクリメントしピクセル数と比較、終了ではない場合は読み込みと書き込みのメモリアドレスをそれぞれインクリメントします。この処理の繰り返しで変換が完了します。
NEONを使うように実装
次にこの処理をNEONを使用して並列化させることを考えます。NEONのレジスタは128bitあるためunsigned char型の8bitのデータであれば一度に16個を保持することができ、NEONの命令セットにもunsigned char型のデータを16個をまとめて1命令で読み込む命令が存在します。さらに今回はRGBの3つの値をまとめて読み込むためインタリーブしつつレジスタに読み込む命令を使用します。これを使用することでメモリ上にRGBRGBRGB…RGBと並ぶデータを3つのレジスタにRRR…R, GGG…G, BBB…Bと1命令で読み込むことができます。よって3要素x16個で合計48バイトを1命令で読み込むことになります。
具体的には”ld3″命令を使用します。ld3命令の動きを図示すると次のようになります。下図では簡単のために8bit x 8個で64bitのデータをイメージしています。
メモリ: R1 G1 B1 R2 G2 B2 R3 G3 B3 R4 G4 B4 R5 G5 B5 R6 G6 B6 R7 G7 B7 R8 G8 B8
メモリアドレスがx1のとき
ld3 {v0 - v2}, [x1], 0 を実行
レジスタv0: R1 R2 R3 R4 R5 R6 R7 R8
レジスタv1: G1 G2 G3 G4 G5 G6 G7 G8
レジスタv2: B1 B2 B3 B4 B5 B6 B7 B8次に読み込んだunsigned char型のデータをfloat型に変換します。データ長を考えると8bitから32bitになります。すでに128bitのレジスタにびっしりと16個の要素を読み込んでいるので16個のまま32bitのデータを保持するためには128bitレジスタが4つ必要になります。
まず128bitのレジスタ1つに8bitデータが16個入っている状態を考えます。これを16bitデータ16個にするには次のような命令を使用します。ushll命令で先頭の8つを16bitに変換し、ushll2で後ろの8個を16bitに変換します。
ushll v1.8h, v20.8b, 0
ushll2 v0.8h, v20.16b, 0次に16bitを32bitに変換します。上と同様の手順で行います。
ushll v3.4s, v1.4h, 0
ushll2 v1.4s, v1.8h, 0次に32bitのデータを浮動小数点のデータに変換します。これはscvtf命令を使用します。
scvtf v1.4s, v1.4sこのようにしてメモリから16ピクセル分をまとめて読み込み、8bitから32bitの浮動小数点に変換することで(32bit/8bit * 3(RGB) ) = 12レジスタ分のデータが用意されます。
この動きを図示すると次のようになります。ここでは0からfまでの8bitの値が16bit, 32bitと順に拡張されていくことを表しています。
レジスタv0: [ 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | a | b | c | d | e | f ], 8bit x 16
ushll v1.8h, v0.8b, 0
ushll2 v2.8h, v0.16b, 0 とすると
レジスタv1: [ 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 ], 16bit x 8
レジスタv2: [ 8 | 9 | a | b | c | d | e | f ], 16bit x 8
ushll v3.4s, v1.4h, 0
ushll2 v4.4s, v1.8h, 0
ushll v5.4s, v2.4h, 0
ushll2 v6.4s, v2.8h, 0 とすると
レジスタv3: [ 0 | 1 | 2 | 3 ], 32bit x 4
レジスタv4: [ 4 | 5 | 6 | 7 ], 32bit x 4
レジスタv5: [ 8 | 9 | a | b ], 32bit x 4
レジスタv6: [ c | d | e | f ], 32bit x 4次は用意したデータに変換用の係数を掛けます。係数は定数なのでリードオンリーなデータとして配置するようにします。赤色に対応する係数は0.299なのでこれのバイナリ値である0x3E991687を設定します。vレジスタとqレジスタが同じ領域を指すことを利用して下記のように行います。
adrp x3, Lred
ldr q7, [x3, #:lo12:.Lred]
.LFE0:
.section .rodata
.align 4
.Lred:
.word 0x3e991687
.word 0x3e991687
.word 0x3e991687
.word 0x3e991687同様に緑や青の係数もロードします。
レジスタに係数をロードしたらそれらとの積をとります。浮動小数点の積をとるためにはfmul命令を使用します。
fmul v3.4s, v3,4s, v7.4s
fmul v1.4s, v1.4s, v7.4s
fmul v2.4s, v2.4s, v7.4s
fmul v0.4s, v0.4s, v7.4s赤色成分の値を求めたら緑と青の値はそこに積をとりつつ加算していきます。積和命令のfmla命令を使用します。緑色の計算の例は次の通りです。
fmla v1.4s, v18.4s, v6.4s
fmla v0.4s, v17.4s, v6.4s
fmla v3.4s, v22.4s, v6.4s
fmla v2.4s, v21.4s, v6.4sこのように処理し赤、緑、青の重み付き和を求めたら32bit浮動小数点を整数値に変換し8bitのデータに変換します。
32bitから16bit, 16bitから8bitのように変換するにはxtnとxtn2命令を使用し、浮動小数点から整数値に変換するにはfcvtzs命令を使用します。
v0, v1, v2, v3レジスタに値が入っている場合の変換例は次の通りです。
fcvtzs v0.4s, v0.4s
fcvtzs v1.4s, v1.4s
fcvtzs v2.4s, v2.4s
fcvtzs v3,4s, v3,4s
xtn v4.4h, v3.4s
xtn2 v4.8h, v1.4s
xtn v1.4h, v2.4s
xtn2 v1.8h, v0.4s
xtn v0.8b, v4.8h
xtn2 v0.16b, v1.8hここまででv0レジスタに変換した8bitのデータが16個入っているのでstr命令でメモリに書き込みます。
str q0, [x3], 16この一連の処理を画像データ全体に適用してRGBからYへの変換は完了です。実際には16ピクセルずつの処理なので16で割った余りが発生する場合は別途処理が必要になります。(NEONを使用しない並列化なしの処理を行うなど)
NEONを使用するように実装したコード
上記の処理をまとめて実装したコードがこちらになります。
.text
.align 2
.p2align 4,,15
.global convertRGBToY
convertRGBToY:
.Lstart:
cbz x2, .Lexit
sub x3, x2, #1
cmp x3, 14
bls .Ljump_init_amari
adrp x3, .Ldata_red
ldr q31, [x3, #:lo12:.Ldata_red]
adrp x3, .Ldata_green
ldr q30, [x3, #:lo12:.Ldata_green]
adrp x3, .Ldata_blue
ldr q29, [x3, #:lo12:.Ldata_blue]
and x5, x2, -16
add x5, x5, x1
mov x4, x0
mov x3, x1
.p2align 3
.Lconv_neon:
ld3 {v20.16b - v22.16b}, [x4], 48
# red
ushll v16.8h, v20.8b, 0
ushll2 v17.8h, v20.16b, 0
ushll v4.4s, v16.4h, 0
ushll2 v5.4s, v16.8h, 0
ushll v6.4s, v17.4h, 0
ushll2 v7.4s, v17.8h, 0
# green
ushll v16.8h, v21.8b, 0
ushll2 v17.8h, v21.16b, 0
ushll v8.4s, v16.4h, 0
ushll2 v9.4s, v16.8h, 0
ushll v10.4s, v17.4h, 0
ushll2 v11.4s, v17.8h, 0
# blue
ushll v16.8h, v22.8b, 0
ushll2 v17.8h, v22.16b, 0
ushll v12.4s, v16.4h, 0
ushll2 v13.4s, v16.8h, 0
ushll v14.4s, v17.4h, 0
ushll2 v15.4s, v17.8h, 0
# convert from integer to float
scvtf v4.4s, v4.4s
scvtf v5.4s, v5.4s
scvtf v6.4s, v6.4s
scvtf v7.4s, v7.4s
scvtf v8.4s, v8.4s
scvtf v9.4s, v9.4s
scvtf v10.4s, v10.4s
scvtf v11.4s, v11.4s
scvtf v12.4s, v12.4s
scvtf v13.4s, v13.4s
scvtf v14.4s, v14.4s
scvtf v15.4s, v15.4s
# multiply and add
fmul v0.4s, v4.4s, v31.4s
fmul v1.4s, v5.4s, v31.4s
fmul v2.4s, v6.4s, v31.4s
fmul v3.4s, v7.4s, v31.4s
fmla v0.4s, v8.4s, v30.4s
fmla v1.4s, v9.4s, v30.4s
fmla v2.4s, v10.4s, v30.4s
fmla v3.4s, v11.4s, v30.4s
fmla v0.4s, v12.4s, v29.4s
fmla v1.4s, v13.4s, v29.4s
fmla v2.4s, v14.4s, v29.4s
fmla v3.4s, v15.4s, v29.4s
# convert from float to integer
fcvtzs v24.4s, v0.4s
fcvtzs v25.4s, v1.4s
fcvtzs v26.4s, v2.4s
fcvtzs v27.4s, v3.4s
xtn v16.4h, v24.4s
xtn2 v16.8h, v25.4s
xtn v17.4h, v26.4s
xtn2 v17.8h, v27.4s
xtn v28.8b, v16.8h
xtn2 v28.16b, v17.8h
str q28, [x3], 16
cmp x3, x5
bne .Lconv_neon
and x6, x2, -16
cmp x2, x6
beq .Lexit
.p2align 3
.Linit_amari:
add x4, x6, x6, lsl 1
add x0, x0, x4
mov w3, 0x45a2
movk w3, 0x3f16, lsl 16
fmov s5, w3
mov w3, 0x1687
movk w3, 0x3e99, lsl 16
fmov s4, w3
mov w3, 0x78d5
movk w3, 0x3de9, lsl 16
fmov s3, w3
.p2align 3
.Lamari:
ldrb w5, [x0, 1]
ldrb w4, [x0, 0]
ldrb w3, [x0, 2]
add x0, x0, 3
ucvtf s2, w5
ucvtf s1, w4
ucvtf s0, w3
fmul s2, s2, s5
fmadd s1, s1, s4, s2
fmadd s0, s0, s3, s1
fcvtzu w3, s0
strb w3, [x1, x6]
add x6, x6, 1
cmp x2, x6
bhi .Lamari
.Lexit:
ret
.p2align 3
.Ljump_init_amari:
mov x6, 0
b .Linit_amari
.Ldata:
.section .rodata
.align 4
.Ldata_red:
.word 0x3e991687
.word 0x3e991687
.word 0x3e991687
.word 0x3e991687
.Ldata_green:
.word 0x3f1645a2
.word 0x3f1645a2
.word 0x3f1645a2
.word 0x3f1645a2
.Ldata_blue:
.word 0x3de978d5
.word 0x3de978d5
.word 0x3de978d5
.word 0x3de978d5解説するほどのことでもありませんが、プログラム中の埋め込みのバイナリは0x 3e99 1687 = 0.299
0x 3f16 45a2 = 0.587
0x 3de9 78d5 = 0.114
を表す数値です。
Linit_amariラベルからの処理で16バイトで割った余りのピクセルを処理しています。ここではNEONを使用しない1バイトごとのアクセスで処理しています。
16で割った余りを求める意味で-16とのANDを求める処理がいくつかあります。これは-16 = 0xff ff ff f0 (32bitの場合)となることを利用して最下位ビットから8bitを0でクリアしています。
また、p2alignはアライメントを指定する書き方で、例えばp2align 3なら2^3=8バイトアライメントされます。
C言語からアセンブリで実装した関数を呼び出すとしっかり動くことが確認できました。試しに適当な画像を入力してみた結果が下の図です。左のカラー画像が入力画像で右の白黒画像が出力結果です。いい感じに白黒になっているのが確認できます。
どれくらい高速化したか
C言語で実装したものとアセンブリでNEONを使用したもので処理速度を比較してみます。白黒に変換する関数を作成し、それを呼び出すことと時間計測するメインルーチンからなるプログラムを作成しました。C言語の最適化は変換の関数のみに適用し、メインルーチンは変換の関数によらず共通となるようにしています。
動作環境はラズパイ3で、64bitカーネルを使用します。測定方法は各100回ずつ行い集計します。時間の単位はマイクロ秒です。一応ナノ秒の精度があるC言語のtimespec_get()を使用して時間取得しています。
まずは比較的小さめの画像サイズで実験しました。入力する画像は28×28ピクセルです。
| 実装方法 | 最適化のレベル | 平均時間 | 中央値 | 標準偏差 |
| C言語 | O2 | 24.309 | 21.770 | 16.164 |
| C言語 | Ofast | 23.052 | 21.589 | 8.697 |
| アセンブラ | – | 4.327 | 4.063 | 1.187 |
次に大きめの画像で実験しました。画像サイズは1620×1080ピクセルです。
| 実装方法 | 最適化のレベル | 平均時間 | 中央値 | 標準偏差 |
| C言語 | O2 | 45,313.038 | 45,308.739 | 62.052 |
| C言語 | Ofast | 45,304.688 | 45,288.323 | 99.174 |
| アセンブラ | – | 11,150.557 | 11,113.495 | 82.912 |
結果を見ると小さい画像、大きい画像ともにしっかりと高速化の成果が出ています。小さい画像では5.3倍、大きい画像では4.1倍の高速化になっています。
C言語で実装した処理がどのようにアセンブリされているかgccに-Sオプションをつけて確認してみましたが1ピクセルずつのアクセスでNEONを使用したコードは生成されていませんでした。
gccのオプションをいろいろと試行錯誤してみてgcc -mcpu=cortex-a53 -Ofast -S func.cなどを試しましたがNEONコードを生成させるところまでできませんでした。
追記: C言語からNEONを使用したコードの生成に成功した
C言語で実装した処理をアセンブリさせてもNEONを使用するコードが生成されずに悩んでいましたがうまくいくようになりました。原因はC言語の実装で、いくつか最適化させるために必要な書き方がありました。
NEONが生成されるようになったC言語での実装はこちらです。
#include <stddef.h>
void convertColorSpaceRGBToY(unsigned char* src, unsigned char* dst, size_t len)
{
for(size_t i = 0; i < len; i++) {
dst[i] =
src[i * 3 + 0] * 0.299f +
src[i * 3 + 1] * 0.587f +
src[i * 3 + 2] * 0.114f;
}
}重要だったのが配列の長さを渡している引数size_t lenで、これをunsigned int lenにしてしまうとNEONを使用した最適化がされませんでした。コンパイル時のオプションは次の通りです。
gcc -mcpu=cortex-a53 -Ofast -S func.c90度の回転処理
回転処理の実装もNEONを使用したものを作成しようと思っていましたが、この投稿が思っていたよりも長くなってしまったので別の機会にしようと思います。
付録
グレースケール化をする関数を呼び出すメインルーチン
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include "/root/lotus_1620_1080.c"
extern void convertRGBToY(unsigned char *src, unsigned char *dst, unsigned int len);
int width = 1620;
int height = 1080;
int main(int argc, char **argv)
{
unsigned char *src;
unsigned char *dst;
src = (unsigned char *)malloc(width * height * 3);
dst = (unsigned char *)malloc(width * height * 1);
memcpy(src, gimp_image.pixel_data, width * height * 3);
struct timespec start, end;
timespec_get(&start, TIME_UTC);
convertRGBToY(src, dst, width * height);
timespec_get(&end, TIME_UTC);
printf("%ld, %ld\n", end.tv_sec - start.tv_sec, end.tv_nsec - start.tv_nsec);
free(src);
free(dst);
return 0;
}入力画像のデータは画像エディタのGIMPで生成したC言語データフォーマットのファイルを利用しています(include “….c”の部分とmemcpy()の部分)。時間計測はtimespec_get()を使用しています。
これとアセンブリ実装をまとめてビルドするMakefileの中身は以下になります。
main:
gcc -c func.s -o func.o
gcc func.o main.c -o main時間を集計するコード
時間計測に使用したコードをメモしておきます。まず100回実行して記録を取ることを次のようにやりました。
for n in $(seq 100); do ./main; done > time_neon.csv100回分の時間ログを集計するコードをpythonで作成しました。numpyを使用しています。
import sys
import numpy as np
data = np.loadtxt(sys.argv[1], delimiter=',')
data[:,0] *= (1 * 1000 * 1000) # uint: second -> micro second
data[:,1] /= 1000 # unit: nano second -> micro second
data = data.sum(axis=1)
m = np.mean(data)
md = np.median(data)
s = np.std(data)
print(f"{m:.3f}, {md:.3f}, {s:.3f}")実行は単純に引数にログファイルを指定します。
python calc.py time_neon.csv
コメント