OpenCVを用いて、画像処理の勉強をします。
細かい理論や1からの実装は行いません。

前回では、入力画像の1画素の値を用いて、出力画像の対応する1画素の値を求めるものでした。

今回は、出力画像の1画素の値を求めるために周辺の複数画素を用いる方法を勉強します。

今回も旭川動物園のペンギンの画像を使用して説明していきます。

import cv2
import matplotlib.pyplot as plt

img = cv2.imread(image_path1)
img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
img = cv2.resize(img[250:2250, 2200:4200], dsize=None, fx=0.1, fy=0.1)
plt.imshow(img, 'gray')

空間フィルタリング

注目画素だけでなく、その周辺の画素も含めた領域内の画素値を用いて計算する処理を空間フィルタリングと呼び、そこで用いられるフィルタを空間フィルタと呼びます。空間フィルタは、線形フィルタと非線形フィルタに分けられます。
線形フィルタは入力画像を$f(i,j)$、出力画像は$g(i,j)$とするとき、次の式で計算できます。
$$
g(i,j) = \sum_{n=-W}^W\sum_{m=-W}^Wf(i+m,j+n)h(m,n)
$$
ここで、$h(m,n)$はフィルタの係数を表す配列であり、フィルタの大きさは$(2W+1)×(2W+1)$となります。
これにあてはまらない処理をともなうフィルタは、非線形フィルタになります。
この後で、いくつかのフィルタの紹介をしていきます。

平滑化

画像処理によって滑らかな濃淡変化を与える処理を平滑化と呼びます。
これは、ノイズなどの軽減に用いられたりします。

平均化

領域内の画素値の平均値を求めるフィルタを平均化フィルタと呼びます。
3×3画素の平均化フィルタは次のように表せます。

    \left[\begin{array}{ccc}
    \frac{1}{9}&\frac{1}{9}&\frac{1}{9}\\
    \frac{1}{9}&\frac{1}{9}&\frac{1}{9}\\
    \frac{1}{9}&\frac{1}{9}&\frac{1}{9}
    \end{array}\right]

OpenCVで平均化を行うには、

2D Convolution (画像のフィルタリング)でフィルタを指定する
ブラー処理

があります。

2D Convolution

cv2.filter2D()でフィルタを指定して畳み込み処理を行います。これは、平均化だけでなく他のフィルタリングで使うことができます。
いくつかのカーネルサイズのフィルタで平均化をしてみます。
ddepthは出力画像のデータ型を決める引数で-1にすることで入力画像と同じにできます。

# フィルタを作成
kernel3 = np.ones((3,3),np.float32)/9
kernel5 = np.ones((5,5),np.float32)/25
kernel7 = np.ones((7,7),np.float32)/49

# 畳み込み処理
dst3 = cv2.filter2D(img, ddepth=-1, kernel=kernel3)
dst5 = cv2.filter2D(img, ddepth=-1, kernel=kernel5)
dst7 = cv2.filter2D(img, ddepth=-1, kernel=kernel7)

fig, ax = plt.subplots(1, 3, figsize=(16, 4))
ax[0].imshow(dst3, 'gray');
ax[0].set_xticks([]);
ax[0].set_yticks([]);
ax[0].set_title('kernel=3');

ax[1].imshow(dst5, 'gray');
ax[1].set_xticks([]);
ax[1].set_yticks([]);
ax[1].set_title('kernel=5');

ax[2].imshow(dst7, 'gray')
ax[2].set_xticks([]);
ax[2].set_yticks([]);
ax[2].set_title('kernel=7');

カーネルサイズが大きくなると画像がボケてくるのがわかります。
また、画像のサイズは入力画像と同じになるように補間されています。
この境界外処理方法はborderTypeという引数で設定できます。

ブラー処理

同様の処理はcv2.blur()でできます。
同じカーネルサイズでブラー処理したものを確認してみます。

# カーネルサイズの設定
n = [3, 5, 7]

fig, ax = plt.subplots(1, 3, figsize=(16, 4))
for i in range(3):
    dst = cv2.blur(img, (n[i], n[i]))
    ax[i].imshow(dst, 'gray')
    ax[i].set_title('kernel='+str(n[i]))
    ax[i].set_xticks([])
    ax[i].set_yticks([])

2D Convolutionと同様の画像を出力することができました。

重み付き平均化

平均値を求めるだけでなく、フィルタ原点に近いほど大きな重みを付ける加重平均フィルタというものもあります。
その重みをガウス分布に近づけたものはガウシアンフィルタと呼びます。
次の式で表されます。
$$
h_{g}(x,y) = \frac{1}{2\pi\sigma^2}\exp{(-\frac{x^2+y^2}{2\sigma^2})}
$$

$\sigma$の値が小さいほど平滑化の効果は小さくなり、大きいほど効果が大きくなります。
たとえば3×3のガウシアンフィルタは次のように表せます。

    \left[\begin{array}{ccc}
    \frac{1}{16}&\frac{2}{16}&\frac{1}{16}\\
    \frac{2}{16}&\frac{4}{16}&\frac{2}{16}\\
    \frac{1}{16}&\frac{2}{16}&\frac{1}{16}
    \end{array}\right]

ここでも3つのカーネルサイズで結果を確認します。cv2.GaussianBlur()を使用します。
$\sigma$は一定としています。

# カーネルサイズ・x方向、y方向のsigmaの設定
k = [3, 5, 7]
sx = [2, 2, 2]
sy = [2, 2, 2]


fig, ax = plt.subplots(1, 3, figsize=(16, 4))
for i in range(3):
    dst = cv2.GaussianBlur(img, (k[i], k[i]), sx[i], sy[i])
    ax[i].imshow(dst, 'gray')
    ax[i].set_title('kernel='+str(k[i]))
    ax[i].set_xticks([])
    ax[i].set_yticks([])

見た目では平均化とは変わらないように見えますが、こちらのほうがノイズの除去には良いらしいです。

特定方向の平滑化

平滑化を特定の方向に限って行うフィルタも考えることができます。
次のように左上から右下への対角成分のみ値をもつようなフィルタを考えてみます。

    \left[\begin{array}{ccc}
    \frac{1}{3}&0&0\\
    0&\frac{1}{3}&0\\
    0&0&\frac{1}{3}
    \end{array}\right]

これはcv2.filter2D()を利用して畳み込み処理を行います。

kernel = np.eye(5)/5
dst = cv2.filter2D(img, -1, kernel)
plt.imshow(dst, 'gray')
plt.xticks([]);
plt.yticks([]);

左上から右下にかけてぶれたような画像となり、左上から右下にかけて平滑化していることがわかります。

エッジ抽出

画像中で明るさが急に変化するエッジ部分を取り出すエッジ抽出のためのフィルタについて説明します。

画像を変えて大きな絵馬の画像を使います。

img2 = cv2.imread(image_path2)
img2 = cv2.cvtColor(img2, cv2.COLOR_BGR2RGB)
img2_rgb = cv2.resize(img2[180:1980, 1000:3600, :], dsize=None, fx=0.2, fy=0.2)
plt.imshow(img2_rgb, 'gray')
plt.show()

img2 = cv2.cvtColor(img2_rgb, cv2.COLOR_RGB2GRAY)
plt.imshow(img2, 'gray')

微分フィルタ

デジタル画像における微分は、注目画素と隣接画素との差分となります。
隣接画素を右側とするか左側とするかで差分値は異なります。
例えば、隣接画素を右側とすれば、

    \left[\begin{array}{ccc}
    0&0&0\\
    0&-1&1\\
    0&0&0
    \end{array}\right]

隣接画素を左側とすれば、

    \left[\begin{array}{ccc}
    0&0&0\\
    -1&1&0\\
    0&0&0
    \end{array}\right]

となります。
さらに、左右の差分値の平均をとって注目画素の微分値とする考え方もあります。
その時の微分フィルタは次のように表せます。

    \left[\begin{array}{ccc}
    0&0&0\\
    -\frac{1}{2}&0&\frac{1}{2}\\
    0&0&0
    \end{array}\right]

同様に縦方向の微分フィルタとして例えば、

    \left[\begin{array}{ccc}
    0&\frac{1}{2}&0\\
    0&0&0\\
    0&-\frac{1}{2}&0
    \end{array}\right]

が考えられます。

各画素における横方向の差分を$\Delta_{x}f(i,j)$、縦方向の差分を$\Delta_{y}f(i,j)$としたとき、$(\Delta_{x}f(i,j), \Delta_{y}f(i,j))$は勾配を表す。
また、勾配の大きさ（エッジ強度）は
$$\sqrt{(\Delta_{x}f(i,j))^2 + \Delta_{y}f(i,j))^2}$$
勾配の方向は、
$$\tan^{-1}\frac{\Delta_{y}f(i,j))}{\Delta_{x}f(i,j))}$$
のように表される。

ここでは、微分フィルタを定義してcv2.filter2D()を使って結果を確認します。

img_tmp = img2.astype('int16')

der_x = np.array([[0,0,0],
                  [-1,0,1],
                  [0,0,0]])/2

der_y = np.array([[0,1,0],
                  [0,0,0],
                  [0,-1,0]])/2

# ddepthで出力画像の型を指定
# 微分はマイナスの値も取るので、符号なしにする
der_x_img = cv2.filter2D(img_tmp, ddepth=cv2.CV_16S, kernel=der_x)
der_y_img = cv2.filter2D(img_tmp, ddepth=cv2.CV_16S, kernel=der_y)

der_int = np.sqrt(der_x_img**2 + der_y_img**2)

fig, ax = plt.subplots(1, 3, figsize=(16, 5))

ax[0].imshow(der_x_img, 'bwr');
ax[0].set_title('横方向');
ax[0].set_xticks([]);
ax[0].set_yticks([]);

ax[1].imshow(der_y_img, 'bwr');
ax[1].set_title('縦方向');
ax[1].set_xticks([]);
ax[1].set_yticks([]);

ax[2].imshow(der_int, 'gray');
ax[2].set_title('勾配の大きさ');
ax[2].set_xticks([]);
ax[2].set_yticks([]);

微分に関してはマイナスになることもあるのでデータ型を変えて表示しています。
縦方向、横方向のエッジを抽出していることが確認できます。

プリューウィットフィルタ、ソーベルフィルタ

プリューウィットフィルタ

微分フィルタでは、エッジ部分だけでなくノイズも強調する傾向にあります。
ノイズを抑えながらエッジを抽出するために、例えば縦方向のエッジの抽出では、
横方向に対して微分してから、それと直交する縦方向に関して平滑化処理をします。
この処理全体は以下のように1つのフィルタを施すことと等価となります。
$\frac{1}{6}$の掛け算を除いて、整数値だけで表されるフィルタをプリューウィットフィルタと呼びます。

    \frac{1}{6}
    \left[\begin{array}{ccc}
    -1&0&1\\
    -1&0&1\\
    -1&0&1
    \end{array}\right] =
    \left[\begin{array}{ccc}
    0&\frac{1}{3}&0\\
    0&\frac{1}{3}&0\\
    0&\frac{1}{3}&0
    \end{array}\right]
    \left[\begin{array}{ccc}
    0&0&0\\
    -\frac{1}{2}&0&\frac{1}{2}\\
    0&0&0
    \end{array}\right]

ここでも、フィルタを定義してcv2.filter2D()で畳み込み処理の結果を確認します。

img_tmp = img2.astype('int16')

kernely = np.array([[-1, 0, 1],
                   [-1, 0, 1],
                   [-1, 0, 1]])/6
kernelx = np.rot90(kernely, -1)

dst1 = cv2.filter2D(img_tmp, ddepth=cv2.CV_16S, kernel=kernely)
dst2 = cv2.filter2D(img_tmp, ddepth=cv2.CV_16S, kernel=kernelx)
dst3 = np.sqrt(dst1**2 + dst2**2)

fig, ax = plt.subplots(1, 3, figsize=(16, 5))
ax[0].imshow(dst1, 'bwr');
ax[0].set_xticks([]);
ax[0].set_yticks([]);
ax[0].set_title('横方向')

ax[1].imshow(dst2, 'bwr');
ax[1].set_xticks([]);
ax[1].set_yticks([]);
ax[1].set_title('縦方向')

ax[2].imshow(dst3, 'gray')
ax[2].set_xticks([]);
ax[2].set_yticks([]);
ax[2].set_title('勾配の大きさ')

ソーベルフィルタ

平滑化の時、中央に重みを付けて平均化する方法もあります。処理全体は以下に示すフィルタと等価となります。
この$\frac{1}{8}$の掛け算を除いた部分のフィルタは、ソーベルフィルタと呼ばれます。

    \frac{1}{8}
    \left[\begin{array}{ccc}
    -1&0&1\\
    -2&0&2\\
    -1&0&1
    \end{array}\right] =
    \left[\begin{array}{ccc}
    0&\frac{1}{4}&0\\
    0&\frac{2}{4}&0\\
    0&\frac{1}{4}&0
    \end{array}\right]
    \left[\begin{array}{ccc}
    0&0&0\\
    -\frac{1}{2}&0&\frac{1}{2}\\
    0&0&0
    \end{array}\right]

まずは、フィルタの定義を行いcv2.filter2()で畳み込み処理をした結果を見てみます。

kernely = np.array([[-1, 0, 1],
                   [-2, 0, 2],
                   [-1, 0, 1]])/8
kernelx = np.rot90(kernely, -1)

dst1 = cv2.filter2D(img_tmp, ddepth=cv2.CV_16S, kernel=kernely)
dst2 = cv2.filter2D(img_tmp, ddepth=cv2.CV_16S, kernel=kernelx)
dst3 = np.sqrt(dst1**2 + dst2**2)

fig, ax = plt.subplots(1, 3, figsize=(16, 5))
ax[0].imshow(dst1, 'bwr');
ax[0].set_xticks([]);
ax[0].set_yticks([]);
ax[0].set_title('横方向')

ax[1].imshow(dst2, 'bwr');
ax[1].set_xticks([]);
ax[1].set_yticks([]);
ax[1].set_title('縦方向')

ax[2].imshow(dst3, 'gray')
ax[2].set_xticks([]);
ax[2].set_yticks([]);
ax[2].set_title('勾配の大きさ')

ソーベルフィルタはcv2.Sobel()が用意されているので、こちらの結果も確認します。

dst1 = cv2.Sobel(img_tmp, ddepth=cv2.CV_16S, dx=1, dy=0, ksize=1)
dst2 = cv2.Sobel(img_tmp, ddepth=cv2.CV_16S, dx=0, dy=1, ksize=1)
dst3 = np.sqrt(dst1**2 + dst2**2)

fig, ax = plt.subplots(1, 3, figsize=(16, 5))
ax[0].imshow(dst1, 'bwr');
ax[0].set_xticks([]);
ax[0].set_yticks([]);
ax[0].set_title('横方向')

ax[1].imshow(dst2, 'bwr');
ax[1].set_xticks([]);
ax[1].set_yticks([]);
ax[1].set_title('縦方向')

ax[2].imshow(dst3, 'gray')
ax[2].set_xticks([]);
ax[2].set_yticks([]);
ax[2].set_title('勾配の大きさ')

微分フィルタと比べてノイズが抑えられ、より滑らかにエッジが抽出されているのがわかります(？)

2次微分とラプラシアン

今までの処理は1次微分に相当するものでした。次に2次微分に対応する処理を説明します。
2次微分は微分を2回繰り返すことであり、デジタル画像では差分を2回繰り返せばよいことがわかります。
注目画素の右と左に半画素ずれた位置の差分値を求めれば、ちょうど注目画素位置の2次微分を求めることになります。
この処理全体は次に示すように1つのフィルタで表すことができて、これを2次微分フィルタと呼びます。

    \left[\begin{array}{ccc}
    0&0&0\\
    1&-2&1\\
    0&0&0
    \end{array}\right] =
    \left[\begin{array}{ccc}
    0&0&0\\
    0&-1&1\\
    0&0&0
    \end{array}\right] - 
    \left[\begin{array}{ccc}
    0&0&0\\
    -1&1&0\\
    0&0&0
    \end{array}\right]

一般に、関数$f(x,y)$のラプラシアンは以下で定義されます。
$$
\frac{\partial^2}{\partial x^2}f(x,y) + \frac{\partial^2}{\partial y^2}f(x,y)
$$
縦方向と横方向の2次微分の結果を足し合わせることでラプラシアンの値を計算できます。
この処理は以下のように1つのフィルタで表すことができて、これをラプラシアンフィルタと呼びます。

    \left[\begin{array}{ccc}
    0&1&0\\
    1&-4&1\\
    0&1&0
    \end{array}\right] =
    \left[\begin{array}{ccc}
    0&0&0\\
    1&-2&1\\
    0&0&0
    \end{array}\right] + 
    \left[\begin{array}{ccc}
    0&1&0\\
    0&-2&0\\
    0&1&0
    \end{array}\right]

ラプラシアンはcv2.Laplacian()で求めることができます。

dst = cv2.Laplacian(img_gray, ddepth=cv2.CV_16S)
plt.imshow(dst,'bwr');
plt.xticks([]);
plt.yticks([]);

ラプラシアンフィルタは、微分を繰り返すためにノイズを強調してしまいます。
そのために、ガウシアンフィルタで平滑化を行った後に、ラプラシアンフィルタを施す方法があります。
この処理は以下のように1つの式で表すことができます。2次元のガウス分布のラプラシアンを計算すると次のようになります。
$$
h_{log}(x,y) = \frac{x^2+y^2-2\sigma^2}{2\pi\sigma^6}\exp(-\frac{x^2+y^2}{2\sigma^2})
$$
この値を係数とするフィルタをLoG(Laplacian of Gaussian)フィルタと呼びます。
まず、LoGフィルタの処理を行う関数を定義します。

def LoG(img, kx, ky, sigmax, sigmay):
    """LoGフィルタ"""
    imgc = img.astype('float32')
    tmp = cv2.GaussianBlur(imgc, (kx, ky) , sigmax, sigmay)
    dst = cv2.Laplacian(tmp, ddepth=cv2.CV_32F)
    return dst

そして、複数の条件で結果を確認してみます。

k = [3, 3, 5]
sx = [2, 8, 6]
sy = [2, 8, 6]

fig, ax = plt.subplots(1, 3, figsize=(16, 10))
for i in range(3):
    dst = LoG(img2, k[i], k[i] , sx[i], sy[i])
    ax[i].imshow(dst, 'bwr')
    ax[i].set_title('kernal='+str(k[i]))
    ax[i].set_xticks([])
    ax[i].set_yticks([])

LoGフィルタは異なる$\sigma$のガウシアンフィルタに差で近似できることが知られています。
このようなフィルタを、DoG(Difference of Gaussian)フィルタと呼びます。
DoGフィルタは以下の式で定義されます。
$$h_g(x,y:\sigma_{1})-h_g(x,y:\sigma_{2})$$

まず、DoGフィルタの処理を行う関数を定義します。

def DoG(img, kx, ky, sigmax1, sigmay1, sigmax2, sigmay2):
    """DoGフィルタ"""
    img = img.astype('float32')
    tmp1 = cv2.GaussianBlur(img, (kx, ky) , sigmax1, sigmay1)
    tmp2 = cv2.GaussianBlur(img, (kx, ky) , sigmax2, sigmay2)
    dst = tmp1 - tmp2
    return dst

そして、複数の条件で結果を確認します。

k = [3, 3, 5]
s1x = [1, 3, 5]
s1y = [1, 3, 5]
s2x = [3, 11, 7]
s2y = [3, 11, 7]

fig, ax = plt.subplots(1, 3, figsize=(16, 10))
for i in range(3):
    dst = DoG(img2, k[i], k[i] , s1x[i], s1y[i], s2x[i], s2y[i])
    ax[i].imshow(dst, 'bwr')
    ax[i].set_title('kernal='+str(k[i]))
    ax[i].set_xticks([])
    ax[i].set_yticks([])

エッジ検出

空間フィルタリングによるエッジ抽出処理の出力は一般に実数の値を持つ画像なので、エッジを特定するにはさらに処理が必要です。
この処理をエッジ検出と呼びます。

LoGフィルタのゼロ交差によるエッジ検出

ラプラシアンフィルタの結果では、エッジをはさみプラスとマイナスの値が対になって現れています。
エッジの位置を求めるには、値がプラスからマイナス、又はマイナスからプラスへ変化する0になる位置（ゼロ交差）を求めればよいことがわかります。
しかし、単純にLoGが0の位置を検出すると、画像の明るさが一定の領域で誤検出を起こしてしまうので、勾配がある一定以下のものを除外する閾値処理を加えることで誤検出を減らすことができます。

実際にエッジの検出を実装して結果を見てみます。

# LoGフィルタで2次微分を求める
dst1 = LoG(img2, 7, 7, 11, 11)
#  -1から1の値の位置でマスク画像を作成(0となる位置の候補)
dst1 = (dst1 >= -1.5)&(dst1 <= 1.5)
dst1 = dst1.astype('uint8')

# 勾配を計算する
tmp = cv2.GaussianBlur(img2, (3, 3), 2)
dir_x = cv2.Sobel(tmp, ddepth=cv2.CV_16S, dx=1, dy=0, ksize=1)
dir_y = cv2.Sobel(tmp, ddepth=cv2.CV_16S, dx=0, dy=1, ksize=1)
m = np.sqrt(dir_x**2 + dir_y**2)

# 勾配の大きな部分のみを対象としたマスク画像の作成
dst2 = m > 20
dst2 = dst2.astype('uint8')

# dst1とdst2の論理和をとることでエッジの検出を行う
dst = cv2.bitwise_and(dst1, dst2)

fig, ax = plt.subplots(1, 4, figsize=(16, 10))
ax[0].imshow(img2, 'gray');
ax[0].set_title('基画像')
ax[0].set_xticks([]);
ax[0].set_yticks([]);

ax[1].imshow(dst1, 'gray');
ax[1].set_title('LoG(値が0付近)')
ax[1].set_xticks([]);
ax[1].set_yticks([]);

ax[2].imshow(m, 'gray');
ax[2].set_title('勾配(閾値以上)')
ax[2].set_xticks([]);
ax[2].set_yticks([]);

ax[3].imshow(dst, 'gray')
ax[3].set_title('エッジ検出')
ax[3].set_xticks([]);
ax[3].set_yticks([]);

ケニーのエッジ検出

得られるエッジが、なるべく連続で、かつ誤検出エッジを含まないように工夫されたものにケニーのエッジ検出という手法があります。次の3段階でエッジを検出しています。

ガウシアンフィルタでノイズを低減する。標準偏差が大きいとノイズ低減効果は高くなるが、エッジ位置が不正確になる。
次に、画素ごとに縦方向と横方向の微分値を求め、勾配の大きさと方向を計算する。
勾配の方向に沿って勾配の大きさが最大の時だけ、その位置をエッジ候補とする。
勾配の方向$\theta(i,j)$に沿った勾配の大きさの補完値$m_{A}$を$m(i,j-1)$と$m(i+1,j-1)$から、
$m_{B}$を$m(i,j+1)$と$m(i-1,j+1)$からから求め、$m(i,j)>m_{A}$かつ$m(i,j)>m_{B}$のときその位置をエッジ候補とする。
1つの閾値の閾値を設けてエッジの判断をするとエッジがとぎれとぎれになることがある。そこで2つの閾値を設けて処理を行う。

上側閾値$Th_{high}$よりも大きいときはエッジと判断する。
下側閾値$Th_{low}$よりも小さいときはエッジではない判断する。
$Th_{high}$と$Th_{low}$の間の時は、エッジとして検出された画素に隣接しているときだけエッジと判断する。

OpenCVにはcv2.Canny()という関数が用意されているのでこれを使って結果を確認します。
Tlowに設定された値より勾配が小さいときはエッジではないと判断します。
Thighに設定された値より勾配が大きいときはエッジであると判断します。
そして、TlowとThighの間の時は、エッジとして検出された画素に隣接しているときだけエッジと判断します。

Tlow = [0, 15, 30]
Thigh = [30, 60, 120]

fig, ax = plt.subplots(3, 3, figsize=(16, 12))
for l in range(3):
    for h in range(3):
        dst = cv2.Canny(img2, Tlow[l], Thigh[h])
        ax[l][h].imshow(dst, 'gray')
        ax[l][h].set_title('l={} h={}'.format(str(Tlow[l]), str(Thigh[h])))
        ax[l][h].set_xticks([])
        ax[l][h].set_yticks([])

鮮鋭化

元の画像の濃淡を残したままエッジを強調する画像の鮮鋭化処理について説明します。

アンシャープマスキング

処理の流れは次のようになっています。

まず入力画像に対して、ある平滑化処理を施しその結果を元画像から引く。
次に、引き算された画像をある定数倍したうえで、入力画像に足し合わせる。

この一連の処理は1つのフィルタと等価であり、以下のように計算されます。
このようなフィルタは鮮鋭化フィルタと呼ばれます。

    1. \quad
    \left[\begin{array}{ccc}
    0&0&0\\
    0&1&0\\
    0&0&0
    \end{array}\right] -
    \left[\begin{array}{ccc}
    \frac{1}{9}&\frac{1}{9}&\frac{1}{9}\\
    \frac{1}{9}&\frac{1}{9}&\frac{1}{9}\\
    \frac{1}{9}&\frac{1}{9}&\frac{1}{9}
    \end{array}\right] =
    \left[\begin{array}{ccc}
    -\frac{1}{9}&-\frac{1}{9}&-\frac{1}{9}\\
    -\frac{1}{9}&\frac{8}{9}&-\frac{1}{9}\\
    -\frac{1}{9}&-\frac{1}{9}&-\frac{1}{9}
    \end{array}\right]

    2. \quad
    \left[\begin{array}{ccc}
    0&0&0\\
    0&1&0\\
    0&0&0
    \end{array}\right] +
    k × \left[\begin{array}{ccc}
    -\frac{1}{9}&-\frac{1}{9}&-\frac{1}{9}\\
    -\frac{1}{9}&\frac{8}{9}&-\frac{1}{9}\\
    -\frac{1}{9}&-\frac{1}{9}&-\frac{1}{9}
    \end{array}\right] =
    \left[\begin{array}{ccc}
    -\frac{k}{9}&-\frac{k}{9}&-\frac{k}{9}\\
    -\frac{k}{9}&1+\frac{8}{9}k&-\frac{k}{9}\\
    -\frac{k}{9}&-\frac{k}{9}&-\frac{k}{9}
    \end{array}\right]

まず、アンシャープマスキングの処理の関数を定義します。
平滑化フィルタとしてはガウシアンフィルタを使用します。

def unsharp_masking(img, kx, ky, sigx, sigy, k):
    img_copy = img.astype('int16').copy()
    img_mean = cv2.GaussianBlur(img_copy, (kx, ky) , sigx, sigy)
    diff_img = img_copy -img_mean
    img_k = diff_img * k
    result = img_copy + img_k
    return result

定義された関数を用いて処理の結果を確認します。

k = [3, 6, 12]

fig, ax = plt.subplots(2, 3, figsize=(16, 8))
for i in range(3):
    dst = unsharp_masking(img2, 3, 3, 2, 2, k[i])
    ax[0][i].imshow(dst, 'gray', vmax=255, vmin=-255)
    ax[0][i].set_title('k='+str(k[i]));
    ax[0][i].set_xticks([]);
    ax[0][i].set_yticks([]);
    
    dst = unsharp_masking(img2_rgb, 3, 3, 2, 2, k[i])
    ax[1][i].imshow(dst, 'gray', vmax=255, vmin=-255)
    ax[1][i].set_title('k='+str(k[i]));
    ax[1][i].set_xticks([]);
    ax[1][i].set_yticks([]);

フィルタを定義して畳み込み処理で求めることもできます。

e = np.array([[0,0,0],
              [0,1,0],
              [0,0,0]])

a = np.array([[-1,-1,-1],
              [-1, 8,-1],
              [-1,-1,-1]]) / 9

f1 = e + 3*a
f2 = e + 6*a
f3 = e + 12*a

dst1 = cv2.filter2D(img2, ddepth=-1, kernel=f1)
dst2 = cv2.filter2D(img2, ddepth=-1, kernel=f2)
dst3 = cv2.filter2D(img2, ddepth=-1, kernel=f3)

fig, ax = plt.subplots(1, 3, figsize=(16, 5))
ax[0].imshow(dst1, 'gray');
ax[0].set_xticks([]);
ax[0].set_yticks([]);
ax[0].set_title('k=3')

ax[1].imshow(dst2, 'gray');
ax[1].set_xticks([]);
ax[1].set_yticks([]);
ax[1].set_title('k=6')

ax[2].imshow(dst3, 'gray')
ax[2].set_xticks([]);
ax[2].set_yticks([]);
ax[2].set_title('k=12')

1つのフィルタで処理できることがわかりました。

エッジを保持した平滑化

局所領域の選択と平均化を行うフィルタ

局所領域のフィルタを用意し、その中の画素値の分散が最小になる領域を選び、
その領域の平均値を出力とすします。
ここでは各局所領域のフィルタを施した時の平均値画像を確認します。

a = np.array([[1,1,0,0,0],
              [1,1,1,0,0],
              [0,1,1,0,0],
              [0,0,0,0,0],
              [0,0,0,0,0]]) / 7.

b = np.array([[0,1,1,1,0],
              [0,1,1,1,0],
              [0,0,1,0,0],
              [0,0,0,0,0],
              [0,0,0,0,0]]) / 7.

c = np.array([[0,0,0,0,0],
              [0,1,1,1,0],
              [0,1,1,1,0],
              [0,1,1,1,0],
              [0,0,0,0,0]]) / 9.

local_filters = [a, b, np.rot90(a, 3),
                 np.rot90(b, 1), c, np.rot90(b, 3),
                 np.flipud(a), np.flipud(b), np.rot90(a, 2)]

img_mini = img2[240:280, 180:220]

fig, ax = plt.subplots(3, 6, figsize=(16, 8))
for i in range(9):
    ax[i//3][i%3].imshow(local_filters[i], 'gray')
    
    dst = cv2.filter2D(img_mini, -1, local_filters[i])
    ax[i//3][i%3+3].imshow(dst, 'gray')
    ax[i//3][i%3+3].set_xticks([])
    ax[i//3][i%3+3].set_yticks([])

k最近隣平均化フィルタ

注目画素の画素値に対し、N画素×N画素の近傍領域中で近い値を持つk個の画素を選び出し、その平均値を出力とします。
$k<=N^2/2$とすることで、注目領域中の画素数の半分以下の画素の平均化をするため、エッジが保存されやすいです。
このようなフィルタをk最近隣平均化フィルタと呼びます。

バイラテラルフィルタ

バイラテラルフィルタとは、ガウシアンフィルタのような注目画素からの距離に応じた重みに加え、注目画素との画素値の差に応じて同様にガウス分布に従う重みを付けて平均化を行うフィルタです。
入力画像を$f(i,j)$、出力画像を$g(i,j)$とするとき、以下の式で表されます。
$$
g(i,j) = \frac{\sum_{n=-W}^{W}\sum_{n=-W}^{W}w(i,j,m,n)f(i+m,j+n)}{\sum_{n=-W}^{W}\sum_{n=-W}^{W}w(i,j,m,n)}
$$
$$
w(i,j,m,n)=\exp{(-\frac{m^2+n^2}{2\sigma_{1}^2})}\exp{(-\frac{(f(i,j)-f(i+m.j+n))^2}{2\sigma_{2}^2})}
$$
ここで、$w(i,j,m,n)$は画素値$f(i+m,j+n)$に対する重みであり、$\sigma_{1}$、$\sigma_{2}$はそれぞれ空間方向、画素値方向の重みを表すガウス分布の標準偏差です。
この重みは注目画素ごとに変化します。

cv2.bilateralFilter()があるのでこれを利用して結果を確認します。
dは注目画素をぼかすために使われる領域、sigmaColorは色についての標準偏差($\sigma_{2}$)、sigmaSpaceは距離についての標準偏差($\sigma_{1}$)です。

n = [10, 20, 30]
m = [10, 20, 30]

fig, ax = plt.subplots(3, 3, figsize=(16, 12))
for i in range(3):
    for j in range(3):
        dst = cv2.bilateralFilter(img2, d=15, sigmaColor=n[i], sigmaSpace=m[j])
        ax[i][j].imshow(dst, 'gray')
        ax[i][j].set_title('sigmaColor={} sigmaSpace={}'.format(str(n[i]), str(m[j])))
        ax[i][j].set_xticks([])
        ax[i][j].set_yticks([])

ノンローカルミーンフィルタ

ノンローカルミーンフィルタとは注目画素の周りの小領域の画素値パターンと、周辺画素のまわりの小領域の画素値パターンとの間の類似度に基づいた重みを用いた平均化です。
重み付き平均を計算する周辺画素の範囲を大きくしても、画像の劣化があまり生じません。
この重みは以下の式で表されます。
$$
w(i,j,m,n)=\exp{(-\frac{\sum_{t=-W'}^{W'}\sum_{s=-W'}^{W'}(f(i+s,j+t)-f(i+m+s,j+n+r))^2}{2\sigma_{3}^2})}
$$
ここで、$/sigma_{3}$は画素値パターンの類似度に応じた重みを表すガウス分布の標準偏差です。

cv2.fastNlMeansDenoising()を使って結果を確認します。

fig, ax = plt.subplots(1, 2, figsize=(12, 6))

dst1 = cv2.fastNlMeansDenoising(img2,10,10,7,21)
ax[0].imshow(dst1, 'gray')
ax[0].set_xticks([]);
ax[0].set_yticks([]);

dst2 = cv2.fastNlMeansDenoisingColored(img2_rgb,None,10,10,7,21)
ax[1].imshow(dst2)
ax[1].set_xticks([]);
ax[1].set_yticks([]);

メディアンフィルタ

平均値ではなく中央値を出力とするフィルタもあり、これをメディアンフィルタと呼びます。
エッジがあまり影響を受けず、スパイク状のノイズ除去に効果が高いという特徴があります。

cv2.medianBlur()を使って確認します。

n = [3, 5, 7]

fig, ax = plt.subplots(1, 3, figsize=(16, 10))
for i in range(3):
    dst = cv2.medianBlur(img2_rgb, n[i])
    ax[i].imshow(dst)
    ax[i].set_title('kernel='+str(n[i]))
    ax[i].set_xticks([])
    ax[i].set_yticks([])

次回

周波数領域におけるフィルタリング

参考文献