Edited at

行列による畳み込みフィルタリング編 ~Python画像処理の再発明家~

More than 1 year has passed since last update.

画像処理ライブラリに頼らず、行列演算だけで畳みこみフィルタリングをするお話。Pythonistaでも可能

基礎編はこちらから


「再発明家」とは

Open CVとかPillowに頼らず、numpyとmatplotlibを使って、様々な画像処理を実際に書いてみる。iOSアプリのPythonistaでも使える組み合わせだ。

import numpy as np

import matplotlib.pyplot as plt

また、画像表示には以下の関数が便利である。(詳しくは基礎編

def img_show(img : np.ndarray, cmap = 'gray', vmin = 0, vmax = 255, interpolation = 'none') -> None:

'''np.arrayを引数とし、画像を表示する。'''

#dtypeをuint8にする
img = np.clip(img,vmin,vmax).astype(np.uint8)

#画像を表示
plt.imshow(img, cmap = cmap, vmin = vmin, vmax = vmax, interpolation = interpolation)
plt.show()
plt.close()


畳み込みによるフィルタリング

畳み込みによるフィルターは空間フィルタリング がわかりやすい。

まず、2次元配列を畳み込みをするための関数を作成する。

def convolve2d(img, kernel):

#部分行列の大きさを計算
sub_shape = tuple(np.subtract(img.shape, kernel.shape) + 1)

#関数名が長いのでいったん省略
strd = np.lib.stride_tricks.as_strided

#部分行列の行列を作成
submatrices = strd(img,kernel.shape + sub_shape,img.strides * 2)

#部分行列とカーネルのアインシュタイン和を計算
convolved_matrix = np.einsum('ij,ijkl->kl', kernel, submatrices)

return convolved_matrix

上記コードはimgの部分行列の行列を使って畳み込みをしている。詳しくはstackoverflow先生を見てほしい。なお、このアルゴリズムでは周辺部分が削れる。

'tiger.jpeg'をトリミングした画像をNTSC 係数による加重平均法でグレースケール化した画像に対して、様々なフィルターを用いた。

img = plt.imread('tiger.jpeg')[1200:1500,1400:1700]

img = (0.298912 * img[...,0] + 0.586611 * img[...,1] + 0.114478 * img[...,2])
img_show(img)

img.png


ローパスフィルター

ローパスフィルターとは、画像の高周波数部分よりも低周波数部分を残しやすいフィルターのことである。

(ハイ・ローは周波数のこと)

要するに、急激な変化をより乏しくすることで、ぼかしをかけている。この処理はノイズ取りのためによく使われる。

img_show(convolve2d(img, kernel))

下表ではconcatenateを用いて元画像と並べている。

手法名
カーネル・元画像・フィルター画像
備考

単純移動平均フィルター

kernel = np.ones((5,5)/25 sma.png

最も単純なローパスフィルター

ガウスフィルター

kernel = gaussian_kernel(5) gaussian.png

正規分布を基にした優秀なローパスフィルター

元の画像では、中央部の薄灰色(カラー画像では黄緑)の部分が点描になっている。しかし、ローパスフィルターをかけると確認が困難になる。

なお、表中で以下のガウス行列を作る関数を用いている。参考

def gaussian_kernel(n : int) -> np.ndarray:

'''(n,n)のガウス行列を作る'''

#[nC0, nC1, ..., nCn]を作成
combs = [1]
for i in range(1,n):
ratio = (n-i)/(i)
combs.append(combs[-1]*ratio)
combs = np.array(combs).reshape(1,n)/(2**(n-1))

#縦ベクトルと横ベクトルの積でガウス行列を作る
result = combs.T.dot(combs)
return result


ハイパスフィルター

ハイパスフィルターとは、画像の低周波数部分よりも高周波数部分を残しやすいフィルターのことである。

これにより、変化に富んだ部分がより強調され、画像がシャープになる。

一般的な鮮鋭化フィルターは$k>0$に対して

1.\hspace{5pt}

\left(
\begin{matrix}
0&-k&0\\
-k&1+4k&-k\\
0&-k&0\\
\end{matrix}
\right)
\\
\mbox{あるいは}
\\
2.\hspace{5pt}

\left(
\begin{matrix}
-k&-k&-k\\
-k&1+8k&-k\\
-k&-k&-k\\
\end{matrix}
\right)

sharp_kernel_1 = lambda k : np.matrix('0,{0},0;{0},{1},{0};0,{0},0'.format(-k,1+4*k))

sharp_kernel_2 = lambda k : np.matrix('{0},{0},{0};{0},{1},{0};{0},{0},{0}'.format(-k,1+8*k))

とかける。参考

このkの値が大きいほど、強いフィルターになる。また、kが同じなら2のほうが強い。

画像は3*3のガウスフィルターをかけてあり、表示には後述する標準化処理を施している。

img_gaussian = convolve2d(img, gaussian(3))

img_show(norm_img(convolve2d(img_gaussian, kernel)))

下表ではconcatenateを用いて元画像と並べている。

下の薄灰色部分(カラー画像だと黄緑色)は本来点描となっており高周波成分が多い。

特に注目してほしい。

カーネル・元画像・フィルター画像

kernel = sharp_kernel_1(1)sharpen_1_1.png

kernel = sharp_kernel_1(10)sharpen_1_10.png

kernel = sharp_kernel_1(1)sharpen_2_1.png

kernel = sharp_kernel_2(10)sharpen_2_10.png

なお、フィルター画像の画素値は0~255に収まらなくなっているため以下の関数で標準化している。

これは$平均\pm n_{std}\times標準偏差$の区間を0~255に対応させている。(実際には$n_{std}=2$としている)

これは比較のための元画像にもしている処理であり、一番上の画像と多少異なる。

def norm_img(img_filtered,n_std):

'''img_filteredを標準化する'''
mean = img_filtered.mean() #平均
std = img_filtered.std() #標準偏差
img_norm = 256*((img_filtered - mean)/std + n_std)/(n_std*2)
return np.clip(img_norm,0,255).astype(np.uint8)


微分

微分フィルターは輪郭抽出に使われる。というのも、隣のピクセルとの差を求めることで、変化の大きい箇所がわかるからだ。

参考:

実際の結果が表示している

フィルターの違いを様々に説明しているをにした。

画像は3*3のガウスフィルターをかけてある。また、表示にはハイパスフィルターの項に記した標準化処理を施した。

img_gaussian = convolve2d(img, gaussian(3))

img_show(norm_img(convolve2d(img_gaussian, kernel)))

内容
カーネル・元画像・フィルター画像
備考
参考リンク

下から上への単純な一次微分フィルター

kernel = np.matrix('0,-1,0;0,1,0;0,0,0')down_up_1.png

ノイズに弱い
説明付き

下から上へのPrewittフィルター

kernel = np.matrix('-1,-1,-1;0,0,0;1,1,1')down_up_prewitt.png

結果画像付き&説明付き

下から上へのSobelフィルター

kernel = np.matrix('-1,-2,-1;0,0,0;1,2,1')down_up_sobel.png

np.matrix('-1,-1,-1;0,0,0;1,1,1')*1.33との差分は微小

結果画像付き&説明付き

ラプラシアンフィルター

kernel = np.matrix('-1,-1,-1;-1,8,-1;-1,-1,-1')laplacian.png

全方向への二次微分$\nabla^2$

結果画像付き&説明付き


その他

内容
カーネル・元画像・フィルター画像
備考
参考リンク

エンボス加工

kernel = np.matrix('-2,-1,0;-1,1,1;0,1,2')emboss.png

Sobelの中央に1
一番下に結果画像