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

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

**[基礎編はこちらから](http://qiita.com/secang0/items/1229212a37d8c9922901)**

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

```py3
import numpy as np
import matplotlib.pyplot as plt
```

また、画像表示には以下の関数が便利である。（詳しくは[基礎編](http://qiita.com/secang0/items/1229212a37d8c9922901)）

```py3
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()
```

#畳み込みによるフィルタリング
畳み込みによるフィルターは[空間フィルタリング](http://www.clg.niigata-u.ac.jp/~lee/jyugyou/img_processing/medical_image_processing_03_press.pdf) がわかりやすい。


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


```py3
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先生](http://stackoverflow.com/questions/43086557/convolve2d-just-by-using-numpy)を見てほしい。なお、このアルゴリズムでは周辺部分が削れる。

'[tiger.jpeg](http://imgur.com/AfaVcBI)'を[トリミング](http://qiita.com/secang0/items/1229212a37d8c9922901#%E5%9B%9B%E8%A7%92%E5%BD%A2%E3%81%AB%E3%83%88%E3%83%AA%E3%83%9F%E3%83%B3%E3%82%B0330%E8%BF%BD%E5%8A%A0)した画像を[NTSC 係数による加重平均法](http://qiita.com/secang0/items/6fba02d7cc3ef34c0c31#ntsc-%E4%BF%82%E6%95%B0%E3%81%AB%E3%82%88%E3%82%8B%E5%8A%A0%E9%87%8D%E5%B9%B3%E5%9D%87%E6%B3%95)でグレースケール化した画像に対して、様々なフィルターを用いた。



```py3
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](https://qiita-image-store.s3.amazonaws.com/0/172464/e7b74d81-5aaf-4847-e291-527370ee53d3.png)


##ローパスフィルター
ローパスフィルターとは、画像の高周波数部分よりも低周波数部分を残しやすいフィルターのことである。
（ハイ・ローは周波数のこと）
要するに、急激な変化をより乏しくすることで、ぼかしをかけている。この処理はノイズ取りのためによく使われる。

```py3
img_show(convolve2d(img, kernel))
```

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

|手法名  |カーネル・元画像・フィルター画像|備考|
|---|---|---|
|単純移動平均フィルター|`kernel = np.ones((5,5)/25`　![sma.png](https://qiita-image-store.s3.amazonaws.com/0/172464/adb7ed7d-75e5-23ad-57a1-03278ad5469e.png)  |最も単純なローパスフィルター|
|ガウスフィルター|`kernel = gaussian_kernel(5)`　![gaussian.png](https://qiita-image-store.s3.amazonaws.com/0/172464/e287dadf-9007-ad50-1d39-3c63135a5602.png)|正規分布を基にした優秀なローパスフィルター|

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

なお、表中で以下のガウス行列を作る関数を用いている。[参考](http://imagingsolution.blog107.fc2.com/blog-entry-165.html)

```py3
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$に対して

```math
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)
```

```py3
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))
```
とかける。[参考](http://seesaawiki.jp/w/mikk_ni3_92/d/%C1%AF%B1%D4%B2%BD%A5%D5%A5%A3%A5%EB%A5%BF)
このkの値が大きいほど、強いフィルターになる。また、kが同じなら2のほうが強い。

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

```py3
img_gaussian = convolve2d(img, gaussian(3))
img_show(norm_img(convolve2d(img_gaussian, kernel)))
```

下表ではconcatenateを用いて元画像と並べている。
下の薄灰色部分（カラー画像だと黄緑色）は本来点描となっており高周波成分が多い。
特に注目してほしい。

|カーネル・元画像・フィルター画像|
|---|
|`kernel = sharp_kernel_1(1)`![sharpen_1_1.png](https://qiita-image-store.s3.amazonaws.com/0/172464/cd6784d6-9130-0e30-5180-adb934b2d7d0.png)|
|`kernel = sharp_kernel_1(10)`![sharpen_1_10.png](https://qiita-image-store.s3.amazonaws.com/0/172464/b12f452d-e780-e762-09b2-052de21f1c6e.png)|
|`kernel = sharp_kernel_1(1)`![sharpen_2_1.png](https://qiita-image-store.s3.amazonaws.com/0/172464/e69676ef-54e8-6951-3993-e31a534572f7.png)|
|`kernel = sharp_kernel_2(10)`![sharpen_2_10.png](https://qiita-image-store.s3.amazonaws.com/0/172464/e8415e91-062c-321e-dfda-9fdfbe4205b2.png)|

なお、フィルター画像の画素値は0~255に収まらなくなっているため以下の関数で標準化している。
これは$平均\pm n_{std}\times標準偏差$の区間を0~255に対応させている。(実際には$n_{std}=2$としている)
これは比較のための元画像にもしている処理であり、一番上の画像と多少異なる。


```py3
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)
```

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

参考:
**[実際の結果が表示している](http://www.mvision.co.jp/WebHelpIM/_RESOURCE/Filter_01.html)**
**[フィルターの違いを様々に説明している](http://qiita.com/shim0mura/items/5d3cbef873f2dd81d82c)**をにした。

画像は3*3のガウスフィルターをかけてある。また、表示には[ハイパスフィルターの項](http://qiita.com/secang0/items/f3a3ff629988dc660d87#%E3%83%8F%E3%82%A4%E3%83%91%E3%82%B9%E3%83%95%E3%82%A3%E3%83%AB%E3%82%BF%E3%83%BC)に記した標準化処理を施した。


```py3
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](https://qiita-image-store.s3.amazonaws.com/0/172464/8ecf79e4-17fc-bd26-47c3-b87a2ed67d8a.png)|ノイズに弱い|[説明付き](http://qiita.com/shim0mura/items/5d3cbef873f2dd81d82c#1%E6%AC%A1%E5%BE%AE%E5%88%86%E3%83%95%E3%82%A3%E3%83%AB%E3%82%BF)|
|下から上へのPrewittフィルター|`kernel = np.matrix('-1,-1,-1;0,0,0;1,1,1')`![down_up_prewitt.png](https://qiita-image-store.s3.amazonaws.com/0/172464/9e41d380-16e2-bb0a-5047-bd3621e3090d.png)||[結果画像付き](http://www.mvision.co.jp/WebHelpIM/_RESOURCE/Filter_Mvc_Gradient.html)&[説明付き](http://qiita.com/shim0mura/items/5d3cbef873f2dd81d82c#prewitt%E3%83%95%E3%82%A3%E3%83%AB%E3%82%BF)|
|下から上へのSobelフィルター|`kernel = np.matrix('-1,-2,-1;0,0,0;1,2,1')`![down_up_sobel.png](https://qiita-image-store.s3.amazonaws.com/0/172464/9171e92b-fdc6-72ff-63da-7ab5793462c8.png)|`np.matrix('-1,-1,-1;0,0,0;1,1,1')*1.33`との差分は微小|[結果画像付き](http://www.mvision.co.jp/WebHelpIM/_RESOURCE/Filter_Mvc_Sobel.html)&[説明付き](http://qiita.com/shim0mura/items/5d3cbef873f2dd81d82c#sobel%E3%83%95%E3%82%A3%E3%83%AB%E3%82%BF)|
|ラプラシアンフィルター|`kernel = np.matrix('-1,-1,-1;-1,8,-1;-1,-1,-1')`![laplacian.png](https://qiita-image-store.s3.amazonaws.com/0/172464/1d54c245-9763-fdfc-dbba-154fc1182634.png)|全方向への二次微分$\nabla^2$|[結果画像付き](http://www.mvision.co.jp/WebHelpIM/_RESOURCE/Filter_Mvc_Laplacian.html)&[説明付き](http://qiita.com/shim0mura/items/5d3cbef873f2dd81d82c#%E3%83%A9%E3%83%97%E3%83%A9%E3%82%B7%E3%82%A2%E3%83%B3%E3%83%95%E3%82%A3%E3%83%AB%E3%82%BF)|


##その他


|内容|カーネル・元画像・フィルター画像|備考|参考リンク|
|---|---|---|---|
|エンボス加工|`kernel = np.matrix('-2,-1,0;-1,1,1;0,1,2')`![emboss.png](https://qiita-image-store.s3.amazonaws.com/0/172464/6528ff77-0e5b-e0fd-a7ef-04c5e156908f.png)|Sobelの中央に1|[一番下に結果画像](https://docs.gimp.org/ja/plug-in-convmatrix.html)|