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

画像処理ライブラリに頼らず、行列演算だけで畳みこみフィルタリングをするお話。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_show(convolve2d(img, kernel))

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

手法名	カーネル・元画像・フィルター画像	備考
単純移動平均フィルター	kernel = np.ones((5,5)/25	最も単純なローパスフィルター
ガウスフィルター	kernel = gaussian_kernel(5)	正規分布を基にした優秀なローパスフィルター

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

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

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)
kernel = sharp_kernel_1(10)
kernel = sharp_kernel_1(1)
kernel = sharp_kernel_2(10)

なお、フィルター画像の画素値は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')	ノイズに弱い	説明付き
下から上へのPrewittフィルター	kernel = np.matrix('-1,-1,-1;0,0,0;1,1,1')		結果画像付き&説明付き
下から上へのSobelフィルター	kernel = np.matrix('-1,-2,-1;0,0,0;1,2,1')	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')	全方向への二次微分$\nabla^2$	結果画像付き&説明付き

##その他

内容	カーネル・元画像・フィルター画像	備考	参考リンク
エンボス加工	kernel = np.matrix('-2,-1,0;-1,1,1;0,1,2')	Sobelの中央に1	一番下に結果画像