python+OpenCVでエッジを保存した平滑化(BilateralFilter, NLMeansFilter)

以下のようなフィルタを使って、画像の注目画素の周りにある画素値の平均をとってやれば、平滑化されたぼやけた画像が取得できる。

\left(
\begin{matrix}
1/9 & 1/9 & 1/9 \\
1/9 & 1/9 & 1/9 \\
1/9 & 1/9 & 1/9
\end{matrix}
\right)

が、これだと画像が全体的にぼやけてしまうため、エッジも等しくぼやけてしまう。エッジは残しておきたいけど、ノイズも低減したいときにはそれ用のフィルタを利用する。

バイラテラルフィルタ(Bilateral filter)

単純に注目画素の周りにある画素値を平均するのではなく、注目画素の近くにあるものをより重視して反映させようというのが重み付き平均化。じゃあその重みの振り分けをどうするかというときに、正規分布に従って振ればいいんじゃないというのがガウシアンフィルタ。
平均0, 分散ρのガウス分布は以下のように表される。

\frac{1}{ \sqrt{2 \pi \sigma}} \exp
\begin{pmatrix}
- 
\frac{x^2}{2\sigma ^ 2}
\end{pmatrix}

これだと1次元になるので2次元に拡張すると

\frac{1}{ 2 \pi \sigma ^2} \exp
\begin{pmatrix}
- 
\frac{x^2 + y ^2}{2\sigma ^ 2}
\end{pmatrix}

ガウシアンフィルタは注目画素と周辺画素の距離についてガウス関数で近似した重みをかけていたが、バイラテラルフィルタはそれに加えてさらに注目画素と周辺画素の画素値の差についてもガウス関数での重み付けをしている。注目画素との画素値の差が小さい（=同じような色合い、明るさ）のであれば重みが大きくなり、注目画素との画素値の差が大きければ重みは小さくなる。入力画像をf(i,j)、出力をg(i,j)とすると

g(i,j) =
 \frac{ \sum_{n=-w}^{w} \sum_{m=-w}^{w} w(i,j,m,n) f(i+m, j+n) }

{\sum_{n=-w}^{w} \sum_{m=-w}^{w} w(i,j,m,n)}

w(i,j,m,n) = \exp 
\begin{pmatrix}
- \frac{m^2 + n^2}{2 \sigma_{1}^2}
\end{pmatrix}

\exp 
\begin{pmatrix}
- \frac{(f(i,j)-f(i+m, j+n))^2}
{2 \sigma_{2}^2}
\end{pmatrix}

というとんでもない式になる。注目画素と周辺画素の距離の重み付けをwに関する式の前半のexpで、注目画素と計算対象画素の画素値の差をwに関する式の後半のexpで表している。平均化のためのフィルタだけど、wは全部を足して１になるようにはされていないため、カーネルの値を全部足して１になるように分母が必要になる。

参考： http://imagingsolution.net/imaging/bilateralfilter/

cv2.bilateralFilter(src, d, sigmaColor, sigmaSpace[, dst[, borderType]]) → dst
http://docs.opencv.org/3.0-last-rst/modules/imgproc/doc/filtering.html?highlight=laplacian#bilateralfilter

• src: 入力画像
• d: 注目画素をぼかすために使われる領域
• sigmaColor: 色についての標準偏差。これが大きいと、画素値の差が大きくても大きな重みが採用される。
• sigmaSpace: 距離についての標準偏差。これが大きいと、画素間の距離が広くても大きな重みが採用される。

import cv2
from matplotlib import pyplot as plt

img = cv2.imread('images.jpg', cv2.IMREAD_COLOR)
bi = cv2.bilateralFilter(img, 15, 20, 20)
bi2 = cv2.bilateralFilter(bi, 15, 20, 20)
bi3 = cv2.bilateralFilter(bi2, 15, 20, 20)

plt.subplot(2,2,1),plt.imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
plt.title("original")
plt.xticks([]),plt.yticks([])
plt.subplot(2,2,2),plt.imshow(cv2.cvtColor(bi, cv2.COLOR_BGR2RGB))
plt.title("bi")
plt.xticks([]),plt.yticks([])
plt.subplot(2,2,3),plt.imshow(cv2.cvtColor(bi2, cv2.COLOR_BGR2RGB))
plt.title("bi2")
plt.xticks([]),plt.yticks([])
plt.subplot(2,2,4),plt.imshow(cv2.cvtColor(bi3, cv2.COLOR_BGR2RGB))
plt.title("bi3")
plt.xticks([]),plt.yticks([])

plt.show()

複数回フィルターかけてみた結果。3回もかけると、ノイズっぽいものは消えているが、疑似輪郭っぽいものも出ている。ちょっと画像が悪かったかも…。 こっちはノイズの殆ど無い画像にフィルタかけてみたもの。回数を増すごとに、疑似輪郭が強調されてイラストっぽくなっている。

ノンローカルミーンフィルタ(Non-local Means Filter)

バイラテラルフィルタは注目画素の画素値と周辺画素の画素値の差に応じた重みをつけたけど、ノンローカルミーンフィルタはテンプレートマッチングのように周辺画素を含めた領域が、注目画素の周辺領域とどれくらい似通っているかによって重みを決定する。

具体的な画像での解説はここが分かりやすい。
http://opencv.jp/opencv2-x-samples/non-local-means-filter

式は

w(i,j,m,n) = \exp 
\begin{pmatrix}
\frac{
\sum_{t=-w}^{w} \sum_{s=-w}^{w} (f(i+s,j+t) -f(i+m+s, j+n+t))^2
}
{}

\end{pmatrix}

を上の式のg(i,j)に当てはめればいい。
注目画素周りの領域と周辺画素周りの領域の類似度を出してるみたいだけど、正直数式の意味はよく分からない…

cv2.fastNlMeansDenoisingColored(src[, dst[, h[, hColor[, templateWindowSize[, searchWindowSize]]]]])
http://docs.opencv.org/3.0-beta/modules/photo/doc/denoising.html

• src: 入力画像（カラー）

• templateWindowSize: 周辺領域のテンプレートサイズ

• searchWindowSize: 重みを探索する領域サイズ

• h: 輝度成分のフィルタの平滑化の度合い、大きいとノイズが減少するが、エッジ部にも影響する

• hColor: 色成分のフィルタの平滑化の度合い、10にしておけば十分

参考: http://ishidate.my.coocan.jp/opencv310_6/opencv310_6.htm

import cv2
from matplotlib import pyplot as plt

img = cv2.imread('images.jpg', cv2.IMREAD_COLOR)

dst = cv2.fastNlMeansDenoisingColored(img,None,10,10,7,21)


plt.subplot(2,1,1),plt.imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
plt.title("original")
plt.xticks([]),plt.yticks([])
plt.subplot(2,1,2),plt.imshow(cv2.cvtColor(dst, cv2.COLOR_BGR2RGB))
plt.title("NLMeans")

一発でもかなりノイズが減らせている。