画素あたり定数時間ガウシアンフィルタのサーベイ
@sakanazensen君が主催する『Computer Vision Advent Calendar 2012』に自分も微力ながら参戦することにしました。テーマはCV分野でもよく用いられている『ガウシアンフィルタ』についてです。特に処理時間がスケールσに依存しない『（画素あたり）定数時間ガウシアンフィルタ』についてまとめます。これはスケールスペースを活用する手法（SIFT特徴点検出や顕著性マップなど）において役立つかと思われます。私はこのテーマについてここ一年取り組んできたのですが、その中でつかんだコツや実装後の感想などについて主観的にですが列挙してみます。参考文献も挙げますので、興味のある方はぜひトライしてみてください。
ちなみに画像処理では2次元フィルタリングについて考える必要がありますが、ガウシアンフィルタは可分(separable)である（N次元ガウシアンはN個の1次元ガウシアンの積に分解できる）ため、特に言及がない場合は1次元を前提に話をしていますのでご注意ください。
古典的手法
ガウシアンフィルタといえば「畳み込み」することが一般的ですが、スケールσに比例して計算量が増える問題を抱えています。これをなんとかしたいということで、主に信号処理分野において古くから多くのアプローチが研究されてきました。最も古典的な解決法は、空間信号の畳み込み演算を、スペクトル信号の乗算に、FFTを通じて変換することです。乗算は定数時間で実施できるので高速に計算できるようになります。しかし実際にはFFT・逆FFTのコストが無視できないオーバーヘッドとなり、σがかなり大きいケースでないと御利益がありません。そんなこんなで、σに依存しない処理時間で動く近似フィルタ*1が様々提案されてきました。
フィルタの重ねがけ (1980s~)
これは中心極限定理から導かれる「どんなフィルタも何度も重ねがけするとガウシアンフィルタの結果に収束する」という理屈を利用した手法です。どんなフィルタでもよいということで、最もシンプルなボックスフィルタの重ねがけが一番使われています。B-splineフィルタや二項フィルタもこれに該当すると思われます。特に文献[03]はBスプライン近似についてよくまとまっていて全体像を把握しやすいかと思います。近年では、ボックスフィルタの幅などを非整数に拡張することで精度をあげた方法[04]や、可分性(separability)を使わずにN次元を一気に処理するような方法[05]が出てきています。
この方法の利点は計算量です。重ねがけが3〜4回で十分な上、ボックスフィルタ自体も簡単な演算で構成されるので、全体の計算コストは少ないです。
しかし一方で、スケールσの制御が面倒であるという問題があります。欲しいσが具体的に決まっている場合、ボックスフィルタの幅を調整するなどして所望のσへと持っていく必要があるわけですが、この計算が少しややこしいです。また常に重ねがけが要求されるために、例えσが小さくても無理やり重ねがけしなければならず、結局直接畳み込んだほうが速かったというケースもあり得ます。そして残念ながら後述の手法に比べると誤差もやや大きめです。
とにかくそれなりにガウシアンフィルタでありさえすれば、σとかノイズとかあまり気にしない！という場合に向いた方法と思われます。
- [01] M. Unser, A. Aldroubi and M. Eden: "B-spline filter signal processing: Part II. efficient design and applications", IEEE Trans. on SP, vol.41, no.2, pp.843-848 (Feb. 1993).
- [02] M. Aubury and W. Luk: "binomial filters", J. of VLSI Signal Processing (1995).
- [03] Y.P. Wang and S.L. Lee: "scale-space derived from B-splines", IEEE Trans. on PAMI, vol.20, no.10 (Oct. 1998).
- [04] P. Gwosdek, S. Grewenig, A. Bruhn and J. Weickert: "theoretical foundations of Gaussian convolution by extended box filtering", Proc. of SSVM, pp.447-458 (2011).
- [05] J.A. Robinson: "efficient Gaussian filtering using cascaded prefix sums", Proc. of IEEE ICIP, pp.117-120 (Oct. 2012).
再帰型フィルタによる近似 (1990s~)
もう一つのよく知られたアプローチが再帰型フィルタによる近似です。再帰型フィルタとはフィードバック構造を持つフィルタのことで、すなわち既に出力した画素値をそれ以降の入力として再利用します。これにより少ない入力画素数（およそ3〜4項で十分）で巨大なσも再現できるようになります。
再帰型フィルタによるガウシアン近似は二つのタイプに分けられます。並列型[11,12,13]と直列型[14,15]です。この手法の元祖が文献[11,12]、それのオフセットノイズに関する改良が[13]になります。より計算効率の少ない直列型の走りが[14]、その一般化・効率化版が[15]になります。文献[11,12,14,15]は今でもよくリファーされている有名どころですが、お互いの関係性を理解するには向かないかもしれません。文献[16,17]は少しサーベイに近い感じで、導入時にはこちらをオススメします。
フィードバック機構による制限として、再帰型フィルタでは着目画素から見て未来の位置にある画素を使えません。そのためガウシアンカーネルを左右半分ずつに分けて2パスのフィルタ処理することでこの問題を回避しています。このとき左右各パートを並列あるいは直列に処理するかで自由度があります。計算量の観点からは直列型のほうが乗算回数も消費メモリ量も少ないですが、並列処理する場合は逆に欠点となる可能性があります。ただ、2次元である画像処理では各ラインを独立に処理できるため、この欠点はあまり気にしなくてもよいと個人的には思います。
Dericheの提案から20年、再帰型フィルタはまさに state-of-the-art の地位にあり、現在でも医用画像処理など様々な場面で用いられ、ライブラリにも搭載されるなどしています*2。が、細かくみると欠点も見つかります。第一が精度の不安定さです。パラメータにもよりますが私のコードでは32<σでガクッと精度が落ちます（他の論文でも同様の現象が見られます）。おそらく浮動小数の精度が追い付かず、フィードバックが原因で誤差伝搬しているように思われます。第二が画像境界付近での扱いの煩雑さです。フィードバック機構ですので十分に処理ステップを踏むまで出力がガウシアンの形状に収束しません。これを回避するには境界外に余分にバッファを設けて空踏み？させるか、あるいは境界条件を設定した上で方程式を解いて端部分用のフィルタ係数を求める必要があります（このネタだけで数本の論文が見つかる）。
性能そのものは決して悪くないのですが、これらの点について実装時にいろいろ気を配るのは結構面倒だったりもします。印象としては「脚は速いが乗りこなすのが難しいじゃじゃ馬」といった感じです。
- [11] R. Deriche: "fast algorithms for low-level vision", IEEE Trans. on PAMI, vol.12, pp.78-87 (1990).
- [12] R. Deriche: "recursively implementing the gaussian and its derivatives", Research Report 1893, INRIA, France (1993).
- [13] G. Farneback and C.F. Westin: "Improving Deriche-style recursive gaussian filters", J. of Mathematical Imaging and Vision, vol.26, no.3, pp.293-299 (Dec. 2006).
- [14] I.T. Young and L.J. van Vliet: "Recursive implementation of the gaussian filter", Signal Processing, no.44, pp.139-151 (1995).
- [15] L.J. van Vliet, I.T. Young and P.W. Verbeek: "Recursive gaussian derivative filters", Proc. of IAPR ICPR, vol.I, pp.509-514 (Aug. 1998).
- [16] S. Tan, J.L. Dale and A. Johnston: "performance of three recursive algorithms for fast space-variant Gaussian filtering", Real-Time Imaging 9, pp.215-228 (2003).
- [17] D. Hale: "Recursive gaussian filters", CWP-546 (2006).
離散コサイン変換による近似 (2010s~)
近年、新たなスタイルのガウシアンフィルタが出てきました。それは離散コサイン変換と参照テーブル(look-up table)を用いてガウシアンカーネルを高速に畳み込むというものです。ガウシアンカーネルは低周波信号のみで構成されるので、対数時間を消費するFFT的な分割統治法をあえて行わず、低周波成分のみを個別に畳み込んでいくというのが高速化の肝になります。最近流行りの表現をすればスペクトルのスパース性を利用しているとも言えるかもしれません。各周波数の畳み込みはバッファを用いることでワンパスで処理できることが分かっており、メモリアクセス回数の減少、そして高速処理の実現につながります。
文献[21]はこの手法の元祖にあたり、ICIP2011の A Best Student Paper Nominee です。文献[22]はそれを主にメモリアクセスコストの観点から改良し2倍高速化したものです（ただし乗算回数は増加）。文献[23]はさらなる精度・速度の改善であり乗算回数増加の問題を解消しています。・・・気づいた方もいらっしゃるかもしれませんが、これは来月開催のIWAIT＠名古屋大学名古屋工業大学で自分が発表するネタで、その宣伝も兼ねた記事だったりもします（爆）。
利点としては速度と精度、そしてそれらの高い安定性です。上で述べたどの手法よりも高速に走り、その処理速度はσ=2の畳み込みに匹敵します。二次元フィルタリングだと畳み込みに対してさらに高速です。また±3σで truncate されたガウシアンの畳み込みよりも高精度という結果も得られています。理論的にも高校数学の範疇で記述できるので、再帰型フィルタでは難しかった境界条件の導出も容易です。私のテストした限り、少なくともPC上においては上記の中で最高性能を発揮できると思われます。
- [21] E. Elboher and M.Werman: "cosine integral images for fast spatial and range filtering", Proc. of IEEE ICIP, pp.89-92 (Sep. 2011).
- [22] K. Sugimoto and S. Kamata: "Fast image filtering by DCT-based kernel decomposition and sequential sum update", In Proc. of IEEE ICIP, pp.125-128 (Oct. 2012).
- [23] K. Sugimoto and S. Kamata: "Fast Gaussian filter based on DCT-V", In Proc. of IWAIT, six pages (Jan. 2013).
むすび
かなり大雑把でしたが、定数時間ガウシアンフィルタについてまとめてみました。ついでに研究の宣伝もしました（笑）。フィルタ重ねがけと再帰型フィルタについては、有名どころの論文はだいたいカバーしてますので、他の論文を読むときにも話の流れを把握しやすくなるかと思います。
個人的な意見ですが、こういった古くからある定数時間フィルタリングの技術がCV分野ではあまり使われておらず、なんだかもったいないなぁーと思ったことが今回の執筆のキッカケでもあります。今後OpenCVにでもこういった定数時間フィルタ関数が用意されれば、きっとさらなるCVの発展につながるんじゃないかなーと妄想してます。画像解像度もスケール枚数ももっとたくさん使ってやろうじゃありませんかっ！（キリッ。
なお自分のアルゴリズムに関するソースコードは、ご希望があれば来年1月ごろ公開できます（ただし研究用途に限るかもしれません）。
