メディアンフィルタをSIMD化してみた

　ネット検索してもメディアンフィルタ（中央値フィルタ）をSIMD（SSE2などのこと）で高速化しているのをあまり見かけないので作ってみました。
　AviUtl用のフィルタです。いつものところにアップしてあります。

　作った理由は、デジカメを最近買って、ごま塩ノイズを除去するのにはメディアンフィルタが一番だろうと思ったからです。

　さてメディアンフィルタは教科書に載っているような基本的なノイズ除去フィルタです。
　仕組みは注目点およびその周辺の8近傍もしくは4近傍をソートし中央値（メディアン）を求め、それを注目点に置き換えるという単純なものです。
　中央値は平均値と似たような値になるので系統としてはぼかし系のノイズ除去になります。平均値との違いは、極端な値（ごま塩ノイズなど）の影響を受けにくいことです。

【サンプル】
　デジカメは露光時間が長くなると白い点がぽつぽつ出てきます。下左の画像は8秒間露光撮影した画像の一部です。下右は分かりやすいように白い点に赤丸を付けました

メディアンフィルタ(ビフォー）

　処理後は下のとおり。

メディアンフィルタ(アフター)

【分岐しないソート】
　ソートアルゴリズムには様々ありますが、条件分岐があると遅くなるので、分岐しないソートを使うことにしました。
　SSE2には2値を比べてmax（もしくはmin）を求める命令があるので、それを使って比較してスワップということします。組み込み関数を使うとこんな感じです。

#define MIN_MAX_SSE2( A, B ) \
{ \
const __m128i max = _mm_max_epi16( A, B ); \
const __m128i min = _mm_min_epi16( B, A ); \
A = min; \
B = max; \
}

　SSE2未対応CPU用にMMX版もあります。

#define MIN_MAX_MMX( A, B ) \
{ \
const __m64 mask = _mm_cmpgt_pi16( A, B ); \
const __m64 max = _mm_or_si64( _mm_andnot_si64( mask, B ), _mm_and_si64( mask, A ) ); \
const __m64 min = _mm_or_si64( _mm_andnot_si64( mask, A ), _mm_and_si64( mask, B ) ); \
A = min; \
B = max; \
}

　さて、ソートするといっても中央値を求めれば充分なので完全にはソートしません。8近傍（注目点を含め9点）の場合、19回の比較で中央値は求まります。19回の順番はパズルゲームのようなものなので、答えは自分で考えてください。バブルソート風に比較した場合は28回ぐらいだと思います。

　[追記:2010/04/17]
　MMX版は以下のようにXORでやった方が速いようです。楓 software: x86 SIMD Technique アーカイブを参考にしました。

#define MIN_MAX_MMX( A, B ) \
{ \
const __m64 t = _mm_and_si64( _mm_xor_si64( A, B ), _mm_cmpgt_pi16( A, B ) ); \
A = _mm_xor_si64( t, A ); \
B = _mm_xor_si64( t, B ); \
}

　SSE2版も同じように試してみましたが、これは_mm_max_epi*を使用した方が速かったです。もしかしたらCPUによるかもしれません。

【SIMD化】
　上を見れば分かりますが、副作用でSIMD化されてます。
　SSE2はshort型(16bit)を8個同時に計算できるので、平均比較回数は19/8 = 2.375回となります。

【マルチスレッド化】
　処理範囲を分割すれば容易にマルチスレッド化できます。しかもAviUtlにはマルチスレッド用の関数が用意されているので簡単です。
　こういうページも参考になるかと思います。
　並列アプリケーションを作ってみよう - SourceForge.JP Magazine

【高速メモリコピー】
　SIMD化や分岐しないソートで高速化しましたが、ここまでやるとメモリーコピーの速度に近くなります。もっと速くしようとすると高速メモリコピーと同じことをする必要があります。
　高速なメモリーコピーについてはアマレココの作者のこの雑記が参考になります。
　アマレココのサイト：過去雑記2009年12月〜2月
　プリフェッチはあまり効果が無いとあります。たしかにCPUにはハードウェアプリフェッチ機能があるので少し先をプリフェッチしただけでは効果はありません。もし効果を得たいのであればずっと先をプリフェッチする（1ライン先をL2にプリフェッチぐらいのことをする）必要があります。

【ベンチマーク】
　getclk_p.aufで測定。デジカメ用に作ったのでテスト素材はサイズが4000x3000の画像です。有効桁数は3桁。

メディアンフィルタのベンチマーク