画像処理における定番のノイズ除去フィルタとして、バイラテラルフィルタというのがあります。その存在だけは認識してはいるものの、その概要は小難しい式で説明されていることが多く、その全容はよくわかってないです。
というわけで、広く使われているであろうこのバイラテラルフィルタに関して、Python使って実装してみました。簡単に、いい加減に考察してみます。
アルゴリズム
式は書かない。なぜならよくわからないから。アカデミックな資料も多く、論文もあります。その概要を分かりやすく説明しているサイトも多く存在しており、説明を読むと分かった気になります。さらにかみ砕いて説明してみます。
平均化
まずこのノイズが除去できる仕組みですが、値として、周囲の適当な範囲（例えば5x5pixel）の範囲で平均値を取ることで、ノイズを除去する。いわゆる移動平均が分かりやすい。
この平均の代わりにガウス関数を使うのが、いわゆるガウシアンフィルタになります。ガウス分布に沿って、フィルタの重みを変化させているだけで、平均化と本質的には変わらない。以降説明の単純化のため、周辺画素を平均化する事でノイズを除去するものとして説明します。
この平均化の致命的な欠点として、画像が持つ本来のエッジもぼかしてしまいます。
色が違うものは仲間外れ
この適当な範囲で平均化する際に、中心画素値と色が違う画素（値が遠い画素）はその平均化の仲間から外す。というやり方で、画像が持つ本来のエッジをぼかさないようにしています。
例えば以下のような、右下への斜めエッジ部ではどうなるかを考えてみます。
左の図のようなエッジでは、明るい側の画素が中心画素にある。その際に、暗い画素を平均化の仲間に入れてしまうと、その平均値は67になる。に対して明るい12画素だけの平均を取ると123と、明るさを残している。今度はその隣の暗い画素が中心に来た右の図では、全員で平均を取ると53になるのに対し、同様に暗い16画素で平均を取ると15となり、やはり暗さを保存できている。
このようにエッジの先鋭性を残しながら、周囲の画素の平均化ができており、ノイズの除去が可能になっている。
値が違うもの（遠いもの）は仲間外れにして平均化することにより、エッジを残しながらノイズを除去することができています。
仲間外れはさすがにやりすぎ
この仲間外れ方式だと、近くも遠くもない色の扱いが難しくなります。仲間に入れるべきか、入れざるべきか判断がつきにくいです。そこでここにも重みを導入します。この近くも遠くもない色は、半分だけ仲間に入れる。近いものには重く、遠いものには軽い重みをかけてあげた後、平均化することで、いい塩梅にしてあげます。
この半分仲間に入れるというのを、中心画素との色の差（画素値の差）に応じて、フィルタ係数に重みをかけることでこれを実現させています。
この重みにもガウス分布を使います。なのでガウス関数が2回出てきます。この辺りが式の理解を放棄させる理由の一つにもなりそうです。
Pythonでバイラテラルフィルタを実装してみる
もちろんOpenCV使えば1行で終わってしまいますが、それじゃあつまらないので、自前で似たような処理を書いてみます。5×5のバイラテラルフィルタもどきです。
コード
import cv2
import numpy as np
import matplotlib.pyplot as plt
import itertools
img = cv2.imread("input.jpg").astype(np.int32) #符号付きdiffを取るのでsinged 32bit
lut = np.arange(0,256,1)
#標準偏差0.1のガウス関数(255で正規化)
lut = (255*np.exp(-(lut/255)**2 / (2*0.1**2))).astype(np.uint8) #cv2.LUT使うのでunsigned 8bit
#plt.figure() #LUT確認用
#plt.plot(lut)
#ガウシアンフィルタ（16倍して使用）
gaus = np.array([[1,4,6,4,1],
[4,16,24,16,4],
[6,24,36,24,6],
[4,16,24,16,4],
[1,4,6,4,1]])*16
center = img[2:-3,2:-3,:] #中心位置（注目座標）
sums = np.zeros_like(center) #合計値保存用
divs = np.zeros_like(center) #剰数保存用
for x, y in itertools.product(range(5), range(5)):
dif = cv2.LUT(abs(center - img[y:y-5,x:x-5,:]).astype(np.uint8),lut).astype(np.int32)
gau = gaus[x][y]*dif//255
sums += gau*img[y:y-5,x:x-5,:]
divs += gau
cv2.imwrite("output.jpg",sums/divs)
この処理の結果が以下です。まぁなんとなく近しい感じの結果が得られています。大きく間違ってはいなそうです。もっとクールに書けるかもしれませんが。
左上がノイズ入り入力、右上が単純なガウシアンフィルタ、左下が今回の、右下が適当なパラメータで行ったOpenCVのバイラテラル。ガウシアンフィルタがボケボケなのに対して、下２つの結果が良好。
解説
まず最初にテーブル類を作っています。中心から正の方向へのガウス関数をあらかじめテーブルとして作っておきます。8bit unsignedに丸めときます。中心画素と周りの画素との差の絶対値を、このテーブルで変換してあげることで、重みを算出します。
ここでガウス関数のテーブルを事前に作っておく理由は、全画素毎にまじめにnp.exp()を呼んでたらどえらい処理時間がかかるので、先に計算結果を入れておくルックアップテーブルを用意しておき、それを引く形で時間節約するためです。
プロットするとこんな感じ。ガウス分布の中心より右側だけのテーブルです。
実際こんな凝った形じゃなく、線形でもいいじゃないかと思ってしまいますが、ガウスでも線形でもLUTにしてしまえば処理時間は同じなので、いったん教科書通りガウスにします。
また畳みこむフィルタ自体は、合計256のこんな感じのよく見る形。
#ガウシアンフィルタ（16倍して使用）
gaus = np.array([[1,4,6,4,1],
[4,16,24,16,4],
[6,24,36,24,6],
[4,16,24,16,4],
[1,4,6,4,1]])*16
ビット落ち防ぐために16倍して使ってます。
for x, y in itertools.product(range(5), range(5)):
dif = cv2.LUT(abs(center - img[y:y-5,x:x-5,:]).astype(np.uint8),lut).astype(np.int32)
gau = gaus[x][y]*dif//255
sums += gau*img[y:y-5,x:x-5,:]
divs += gau
このfor文で、25画素の差の絶対値を取って、LUT通して、キャストして、重みを乗じて、フィルタ係数画素に掛けて、といった一連の処理を行ってます。
本当は中心画素（x=2,y=2）では差を取る必要ない（必ず0）なので、そこでif文入れた方が速いとは思いましたが、可読性が悪くなるのでここでは無視しました。
ここで、差に応じたガウス関数の値を求めるために、LUTを使っています。OpenCVのLUTを使いたかったので、8bitのガウステーブルにしてます。だいぶビット落ちてます。出力floatのままLUT使いたいです。どうやるんだかわかりませんでした。
加えて、フィルタの係数に重みを掛けたうえでキャストしちゃってるので、ビット欠損が起こっています。一応それを回避するためにフィルタを16倍してますが、足りてるかな？あまり考察してないです。
ほんとは画素値を乗じてから重み掛けた方が欠損は少ないです。あぁなんてことをと思われるかもしれませんが、最後係数の総和で除算するときに値が回り込まないようにするためにはしょうがないかと思ってます。画素毎に剰数が変わるのはめんどくさいです。
さらに、256で割りたい（8bitシフトしたい）ところですが、残念ながらそうもなってません。全画素割り算になってしまいます。遅そう。
周囲の2画素も捨ててます。右端に関しては3画素捨ててます。割り切ってます。
最後にフィルタの係数和で除算してファイルに保存しています。
行数としてはこれだけでできています。型変換が随所にあるのが気持ち悪いです。まぁこんなもんか。
まとめ
今回、自分が分かりやすいと思ったサイトから、さらにかみ砕いて、仲間外れ付き平均という概念で、適当に解説を入れてみました。併せてPythonで実装もしてみました。
一応理解は大きくは間違ってはいなそうです。式は理解してませんが…。こんなもんだと思います。間違ってたら教えてください。
やはり処理は遅いですね。ちまちまフィルタ係数に重みを掛けたり、除算する値が画素毎に動的に変動したりと、遅そうだし、ハード処理には向かなそう。でも並列化には向きそう。高速化の工夫も色々されているようなので、実際は色々工夫してるんだろうなぁ。
追記
少し真面目なコードも書いたので、こっちを見てもらえると、もう少し使えるものになっていると思います。
コメント