はじめに
入門記事の前に、開発者がSIMD命令を使用して画像処理の最適化を最適化する場合に直面しなければならない問題のリストを挙げました。 さて、上記の問題をどのように解決できるかを具体例で示す時が来ました。 最初の例でどのアルゴリズムを選択するかを長い間考えていたため、中央値フィルタリングに焦点を合わせることにしました。 中央値フィルタリングは、暗いシーンでデジタルカメラに必然的に現れるノイズを抑制する効果的な方法です。 このアルゴリズムは非常にリソースを消費します。たとえば、3x3メディアンフィルターでグレー画像を処理する場合、画像ポイントごとに約50回の操作が必要です。 しかし、同時に、彼は8ビットの数値でのみ動作し、動作するのに比較的少ない入力しか必要としません。 これらの状況により、アルゴリズムはSIMD最適化のために十分にシンプルになり、同時に非常に大きな加速を得ることが可能になります。
参考のために、アルゴリズムの本質を思い出します。
- 元の画像の各ポイントに対して、特定の近傍（この場合は3x3）が取得されます。
- この近傍のポイントは、明るさの増加によってソートされます。
- ソートされた近傍の中間点（3x3フィルターでは5番目）が最終画像に記録されます。
最初の例では、タスクを簡素化し、境界効果を考慮せずにグレー画像（1ポイントあたり8ビット）の場合を検討することにしました。このため、わずかに大きな画像の中央部分をフィルタリングに使用しました。 もちろん、実際の問題では、画像の境界点を特別なアルゴリズムを使用して計算する必要があり、これらの特別なアルゴリズムの実装がコードの大部分を占める場合がありますが、これらの問題についてはもう一度説明します。 さらに、同様の理由で、イメージの幅は32の倍数であると想定されます。
アルゴリズムのスカラーバージョン
最初のスカラーバージョン
この最初の実装では、問題は正面から解決されました。元の画像の各点の近傍は補助配列にコピーされ、std :: sort（）関数を使用してソートされました。
inline void Load(const uint8_t * src, int a[3]) { a[0] = src[-1]; a[1] = src[0]; a[2] = src[1]; } inline void Load(const uint8_t * src, size_t stride, int a[9]) { Load(src - stride, a + 0); Load(src, a + 3); Load(src + stride, a + 6); } inline void Sort(int a[9]) { std::sort(a, a + 9); } void MedianFilter(const uint8_t * src, size_t srcStride, size_t width, size_t height, uint8_t * dst, size_t dstStride) { int a[9]; for(size_t y = 0; y < height; ++y) { for(size_t x = 0; x < width; ++x) { Load(src + x, srcStride, a); Sort(a); dst[x] = (uint8_t)a[4]; } src += srcStride; dst += dstStride; } }
2番目のスカラーバージョン
ただし、前の方法ではすぐにボトルネックが発生します。これが標準のソート機能です。 これは、任意のサイズの配列をソートするためのもので、多くの内部チェックを実行するため、配列が小さくサイズが固定されているこの問題の最適な解決策ではない可能性が高くなります。 さらに、配列を完全に並べ替える必要はまったくありません。必要な値が4番目の要素にあれば十分です。 したがって、「ネットワーク」ソート方法（バイトニックソート方法とも呼ばれます）を適用できます。
inline void Sort(int & a, int & b) // { if(a > b) { int t = a; a = b; b = t; } } inline void Sort(int a[9]) // { Sort(a[1], a[2]); Sort(a[4], a[5]); Sort(a[7], a[8]); Sort(a[0], a[1]); Sort(a[3], a[4]); Sort(a[6], a[7]); Sort(a[1], a[2]); Sort(a[4], a[5]); Sort(a[7], a[8]); Sort(a[0], a[3]); Sort(a[5], a[8]); Sort(a[4], a[7]); Sort(a[3], a[6]); Sort(a[1], a[4]); Sort(a[2], a[5]); Sort(a[4], a[7]); Sort(a[4], a[2]); Sort(a[6], a[4]); Sort(a[4], a[2]); }
3番目のスカラーバージョン
最初の方法はゼロ方法よりもはるかに高速であることが判明しましたが（以下のテスト結果を参照）、まだ1つのボトルネックがあります-この方法には多くの条件付き遷移があります。 ご存じのように、最新のプロセッサはパイプラインアーキテクチャを備えており、1クロックサイクルで複数の命令を異常に実行する可能性があるため、あまり好まれていません。 最初と2番目の両方について、条件付き遷移は非常に望ましくありません。これは、プログラムのどのブランチを次に実行するかを見つけるために、プロセッサが停止して計算結果を待つ必要があるためです。 幸いなことに、2つの8ビット値のソートは、条件分岐なしで実装できます。
inline void Sort(int & a, int & b) { int d = a - b; int m = ~(d >> 8); b += d&m; a -= d&m; }
これらの予備的なスカラーコードの最適化の結果は、メディアンフィルターの次のバージョンになります。
inline void Load(const uint8_t * src, int a[3]) { a[0] = src[-1]; a[1] = src[0]; a[2] = src[1]; } inline void Load(const uint8_t * src, size_t stride, int a[9]) { Load(src - stride, a + 0); Load(src, a + 3); Load(src + stride, a + 6); } inline void Sort(int & a, int & b) { int d = a - b; int m = ~(d >> 8); b += d&m; a -= d&m; } inline void Sort(int a[9]) // { Sort(a[1], a[2]); Sort(a[4], a[5]); Sort(a[7], a[8]); Sort(a[0], a[1]); Sort(a[3], a[4]); Sort(a[6], a[7]); Sort(a[1], a[2]); Sort(a[4], a[5]); Sort(a[7], a[8]); Sort(a[0], a[3]); Sort(a[5], a[8]); Sort(a[4], a[7]); Sort(a[3], a[6]); Sort(a[1], a[4]); Sort(a[2], a[5]); Sort(a[4], a[7]); Sort(a[4], a[2]); Sort(a[6], a[4]); Sort(a[4], a[2]); } void MedianFilter(const uint8_t * src, size_t srcStride, size_t width, size_t height, uint8_t * dst, size_t dstStride) { int a[9]; for(size_t y = 0; y < height; ++y) { for(size_t x = 0; x < width; ++x) { Load(src + x, srcStride, a); Sort(a); dst[x] = (uint8_t)a[4]; } src += srcStride; dst += dstStride; } }
このバージョンは、アルゴリズムの初期バージョンよりも大幅に（約5〜6倍）高速です。 そして、それに基づいてアルゴリズムのSIMD最適化を実行し、アルゴリズムのスカラーバージョンとベクトルバージョンの速度を比較します。
アルゴリズムのSIMDバージョン
SIMDを使用したメディアンフィルタリングアルゴリズムの最適化を説明するために、現在は時代遅れであり、より歴史的な関心を集めているMMX拡張を欠いた2組の命令SSE2およびAVX2を使用します。 。 最新のC ++コンパイラのほとんどは、組み込み関数（すべての種類のプロセッサ拡張を使用できる組み込み関数）をサポートしています。 コンパイラ組み込み関数を使用したプログラミングは、実際には純粋なCでのプログラミングと違いはありません。組み込み関数は、コンパイラによってプロセッサ命令に直接変換されますが、プロセッサレジスタを直接操作することはプログラマには見えません。 ほとんどの場合、組み込み関数を使用するプログラムは、アセンブラーで記述されたプログラムよりも速度が劣りません。
SSE2バージョン
SSE2整数コマンドは、ヘッダーファイル<emmintrin.h>で定義されています。 基本タイプは__m128i-128ビットベクトルで、コンテキストに応じて、2 64ビット、4 32ビット、8 x 16ビット、16 x 8ビットの符号付きまたは符号なしの数値のセットとして解釈できます。 ご覧のとおり、これらはベクトル算術演算だけでなく、ベクトル論理演算、およびデータのロードとアンロードのベクトル演算もサポートしています。 以下は、SSE2命令を使用したメディアンフィルターの最適化です。 コードは、私にとってはかなり視覚的です。
#include < emmintrin.h > inline void Load(const uint8_t * src, __m128i a[3]) { a[0] = _mm_loadu_si128((__m128i*)(src - 1)); // 128 16 a[1] = _mm_loadu_si128((__m128i*)(src)); a[2] = _mm_loadu_si128((__m128i*)(src + 1)); } inline void Load(const uint8_t * src, size_t stride, __m128i a[9]) { Load(src - stride, a + 0); Load(src, a + 3); Load(src + stride, a + 6); } inline void Sort(__m128i& a, __m128i& b) { __m128i t = a; a = _mm_min_epu8(t, b); // 2- 8 16 b = _mm_max_epu8(t, b); // 2- 8 16 } inline void Sort(__m128i a[9]) // { Sort(a[1], a[2]); Sort(a[4], a[5]); Sort(a[7], a[8]); Sort(a[0], a[1]); Sort(a[3], a[4]); Sort(a[6], a[7]); Sort(a[1], a[2]); Sort(a[4], a[5]); Sort(a[7], a[8]); Sort(a[0], a[3]); Sort(a[5], a[8]); Sort(a[4], a[7]); Sort(a[3], a[6]); Sort(a[1], a[4]); Sort(a[2], a[5]); Sort(a[4], a[7]); Sort(a[4], a[2]); Sort(a[6], a[4]); Sort(a[4], a[2]); } void MedianFilter(const uint8_t * src, size_t srcStride, size_t width, size_t height, uint8_t * dst, size_t dstStride) { __m128i a[9]; for(size_t y = 0; y < height; ++y) { for(size_t x = 0; x < width; x += sizeof(__m128i)) { Load(src + x, srcStride, a); Sort(a); _mm_storeu_si128((__m128i*)(dst + x), a[4]); // 128 16 } src += srcStride; dst += dstStride; } }
AVX2バージョン
整数のAVX2コマンドは、ヘッダーファイル<immintrin.h>で定義されています。 基本型は__m256i-256ビットのベクトルで、コンテキストに応じて、4つの64ビット、8ビット32、16ビット16、32ビット8ビットの符号付きまたは符号なしの数値のセットとして解釈できます。 AVX2命令セットは主にSSE2命令セットを繰り返しますが（ベクトルの幅が2倍になると）、SSE2には類似物がないかなり多くの命令も含まれます。 以下は、AVX2命令を使用したメディアンフィルターの最適化です。
#include < immintrin.h > inline void Load(const uint8_t * src, __m256i a[3]) { a[0] = _mm256_loadu_si256((__m128i*)(src - 1)); // 256 32 a[1] = _mm256_loadu_si256((__m128i*)(src)); a[2] = _mm256_loadu_si256((__m128i*)(src + 1)); } inline void Load(const uint8_t * src, size_t stride, __m256i a[9]) { Load(src - stride, a + 0); Load(src, a + 3); Load(src + stride, a + 6); } inline void Sort(__m256i& a, __m256i& b) { __m256i t = a; a = _mm256_min_epu8(t, b); // 2- 8 32 b = _mm256_max_epu8(t, b); // 2- 8 32 } inline void Sort(__m256i a[9]) // { Sort(a[1], a[2]); Sort(a[4], a[5]); Sort(a[7], a[8]); Sort(a[0], a[1]); Sort(a[3], a[4]); Sort(a[6], a[7]); Sort(a[1], a[2]); Sort(a[4], a[5]); Sort(a[7], a[8]); Sort(a[0], a[3]); Sort(a[5], a[8]); Sort(a[4], a[7]); Sort(a[3], a[6]); Sort(a[1], a[4]); Sort(a[2], a[5]); Sort(a[4], a[7]); Sort(a[4], a[2]); Sort(a[6], a[4]); Sort(a[4], a[2]); } void MedianFilter(const uint8_t * src, size_t srcStride, size_t width, size_t height, uint8_t * dst, size_t dstStride) { __m256i a[9]; for(size_t y = 0; y < height; ++y) { for(size_t x = 0; x < width; x += sizeof(__m256i)) { Load(src + x, srcStride, a); Sort(a); _mm256_storeu_si256((__m256i*)(dst + x), a[4]); // 256 32 } src += srcStride; dst += dstStride; } }
注：最適化から最大のアクセラレーションを得るには、AVX2命令を含むコードを次のコンパイラオプションでコンパイルする必要があります（/ arch：Visual Studio 2012のAVX、GCCの-march = core-avx2）。 さらに、この場合にプロセッサの動作モードを切り替えると全体的なパフォーマンスに悪影響が及ぶため、コードにAVX2命令とSSE2命令が交互に含まれないことが非常に望ましいです。 上記に基づいて、AVX2およびSSE2バージョンのアルゴリズムを異なる「.cpp」ファイルに配置することをお勧めします。
テスト中
テストは2 MBのグレー画像（1920 x 1080）で実行されました。 時間測定の精度を向上させるために、テストを数回実行しました。 実行時間は、テストの合計実行時間を実行回数で割ることによって得られました。 実行数は、合計実行時間が少なくとも1秒になるように選択されました。これにより、アルゴリズムの実行時間を測定するときに2つの精度の兆候が得られるはずです。 アルゴリズムは、Microsoft Visual Studio 2012上の64ビットWindows 7向けに最大限に最適化されてコンパイルされ、iCore-7 4770プロセッサー（3.4 GHz）で実行されました。
| 実行時間、ミリ秒 | 相対加速度 | ||||
| 最高のスカラー | SSE2 | AVX2 | 最高のスカラー/ SSE2 | 最高のスカラー/ AVX2 | SSE2 / AVX2 |
| 24.814 | 0.565 | 0.424 | 43.920 | 58.566 | 1.333 |
おわりに
テスト結果からわかるように、最新のプロセッサでのSIMD命令の使用による中央値フィルタリングアルゴリズムの加速は、40〜60倍に達する可能性があります。 これは最適化を気にするのに十分な量のようです（もちろん、プログラムで実行速度が重要でない限り）。 この出版物では、使用するコードを可能な限り簡素化しようとしました。
そのため、データのアライメント、境界点の処理などに関連する問題は、範囲を超えて残りました。
これらの質問を次の記事で開示するようにします。 読者がこの機能を実装するバトルコードがどのように見えるかに興味がある場合は、 ここで見つけることができます 。