32
Help us understand the problem. What are the problem?

More than 5 years have passed since last update.

posted at

updated at

行列による畳み込みフィルタリング編 ~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.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 一番下に結果画像
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Sign upLogin
32
Help us understand the problem. What are the problem?