Python3 & OpenCV で画像処理を学ぶ[6] 〜 numpy で畳み込み演算と平滑化フィルタの実装
1. はじめに
今回は、空間フィルタリングの基礎概念に触れて、いくつかの平滑化フィルタを実装します。
理解を深めるため, まず畳み込み演算から自力でnumpyで実装し, 次にOpenCVでのやり方と比較するという過程を踏みます。
2. 基本概念
2.1. 空間フィルタリングとは
これまで触れてきたヒストグラム, トーンカーブ, 描画モードなどは, あくまで「ピクセルごと」のお話でした。つまり, 画像データでは本来, 「どのピクセルがどこにあるか」「隣接するピクセルは何か」などの位置に関わる情報も重要であるはずなのですが, そういったことは一旦無視した話だったというわけです。
そこで, 画像に対して「隣接するピクセルを考慮した変換」をしたくなります。そのような変換を「空間フィルタリング」と呼びます。 具体例としては, 以下のようなものがあります。
- ブラー（画像処理では「平滑化」の一種. 今回扱います）
- メディアン
- 輪郭検出
- シャープネス
2.2. 畳み込み演算とは
そこで行われるのが, 畳み込み Convolution と呼ばれる演算です。これはどのような演算でしょうか。
まずは一次元配列で考えてみます。まず, 入力の配列の他に, カーネルと呼ばれるフィルタ(配列)を以下のような位置に用意します。
次に, カーネルの要素と, 同じ位置にある入力の要素とをそれぞれ掛けて, それらを足し合わせたものを, カーネルの中心の直下の出力の値とします。
( 3×1 + 1×2 + 2×1 = 7 )
そして, カーネルをひとつ隣に動かして, 再び先ほどの計算を行います。これを全マスで行うことで, 出力が得られます。
当然ながら, カーネルが偶数サイズだと真ん中が存在しないので, カーネルは奇数サイズである必要があります。
また, これでは入力の一番端にあたるマスの出力が計算できない気がしてきますが, その場合は以下のような方法で対処します。
- 必要なマスを 0 で埋める(ゼロパディング)
- 最端にある値を繰り返す
- 画面端を境に対称となるように, 反転しながら繰り返す
以上をまとめると, 畳み込みとは以下のような演算ということになります。
画像については, これをそのまま二次元配列に拡張すればよいです。以下のようなイメージですね。
これからまた具体的に見ていきますが, この畳み込み演算で用いるカーネルの値によって, ブラーがかかったりエッジを抽出したりなど, 様々な効果が生まれるということです。
また, 以上のような「カーネルの値と画素値 の積和」という計算によってのみ構成されるフィルタを線形フィルタと呼び, そうでない式を含むフィルタを非線形フィルタと呼ぶようです.
3. 平滑化の実装
import numpy as np from numpy.lib.stride_tricks import as_strided from matplotlib import pyplot as plt from mpl_toolkits.mplot3d import axes3d import matplotlib.gridspec as gridspec import seaborn as sns import cv2 %config InlineBackend.figure_format = 'retina'
3.1. 平均化フィルタ
平均化フィルタは, あるピクセルの周囲の画素値の平均値を出力とするフィルタです。
例えば, 以下のようなカーネルを用います。
kernel = np.ones([3,3]) / (3*3) print(kernel)
[[ 0.11111111 0.11111111 0.11111111]
[ 0.11111111 0.11111111 0.11111111]
[ 0.11111111 0.11111111 0.11111111]]
要するに, 和が 1 になるように, 各マスに均等な値を用います。 5×5なら以下の通りです。
N = 5 # should be odd kernel = np.ones([N,N]) / (N*N) print(kernel)
[[ 0.04 0.04 0.04 0.04 0.04]
[ 0.04 0.04 0.04 0.04 0.04]
[ 0.04 0.04 0.04 0.04 0.04]
[ 0.04 0.04 0.04 0.04 0.04]
[ 0.04 0.04 0.04 0.04 0.04]]
平均とは, 以下のようなものだったことを思い出すと,
$ \dfrac{a_1 + a_2 + … + a_N}{N} $
このようなカーネルによる畳み込み演算の出力は, 周囲の画素値の平均値となっていることがわかります。
3.1.1 python/numpy による実装
さっそく, 畳み込み演算による平均化フィルタを実装してみたいと思います。まずは OpenCV に頼らず, numpy で畳み込み演算を実装してみることにします。
そのために, 先ほど from numpy.lib.stride_tricks import as_strided として読み込んでいた numpy.lib.stride_tricks.as_strided() という関数を用いると便利です。
これは, ある numpy.ndarray から, 要素を抜き出した view を作ることができる関数です。 view とはつまり, 新たな配列を実体として作るのではなく, 元の array が確保したメモリ上の値を参照する配列です。 view の要素の値を書き換えると, 元の array の値も変更されます。
また, numpy.ndarray は, 本来 C の配列のようにメモリ上で連続して確保されています。そこで, この関数を用いることで, 一定 byte 間隔(stride)で値を抜き出してくる, といったことができます。
実装の基本方針としては以下の通りです。
例えば, (1280,720) のimg配列に対する (3,3) のカーネルの畳み込み演算では, as_strided を用いて (1280, 720, 3, 3) の配列を作ります。 (1280,720) の配列が, 畳み込みに備えた (3,3) 配列を要素として持つイメージです。
さて, そうして実装した畳み込み演算が以下になります。
まずはいつも通り, 画像表示用の関数。
def img_show2(img1,img2,titles=['ref_img','filtered']): fig = plt.figure(figsize=(20,7)) ax1 = fig.add_subplot(121) ax2 = fig.add_subplot(122) ax1.set_title(titles[0]), ax2.set_title(titles[1]) ax1.set_xticks([]), ax1.set_yticks([]) ax2.set_xticks([]), ax2.set_yticks([]) ax1.imshow(img1) ax2.imshow(img2) plt.show()
次に, imgとkernelを受け取って畳み込みをした配列を返す関数です。imgが3次元配列だとややこしくなるので, グレースケールもしくはRGB各チャンネルごとに利用することを想定しています。
def convolve2d_1ch(img, kernel, padding='0'): # エラー処理 if len(img.shape) != 2: print('img should be 2d-array') return if kernel.shape[0]%2 == 0: print('Kernel size shold be odd') return edge = int(kernel.shape[0]/2) # 境界をどうするか if padding=='edge': img = np.pad(img, [edge,edge], mode='edge') elif padding=='reflect': img = np.pad(img, [edge,edge], mode='reflect') else: img = np.pad(img, [edge,edge], mode='constant', constant_values=0) sub_shape = tuple(np.subtract(img.shape, kernel.shape) + 1) # (rows, cols, N, N) conv_shape = sub_shape + kernel.shape # A.strides は (Aのrow方向の一要素ごとのbytes数, col方向のbytes数) のタプル # (shape[0]*bytes, bytes, shape[0]*bytes, bytes) strides = img.strides + img.strides submatrices = as_strided(img, conv_shape, strides) # (i,j)配列と(k,l,i,j)配列で, (i,j)要素同士の積和を要素とした(k,l)配列を返す convolved_mat = np.einsum('ij, klij->kl', kernel, submatrices) return convolved_mat
これらを29×29のカーネルを用いて実際に実行すると, 以下のようになります。
img = cv2.cvtColor(cv2.imread('img_src/wizdomiot.png',1), cv2.COLOR_BGR2RGB)
N = 29 kernel = np.ones([N,N]) / (N*N) mean_convolved = np.zeros(img.shape) for (i,ch) in enumerate(['R','G','B']): mean_convolved[:,:,i] = convolve2d_1ch(img[:,:,i], kernel, padding='edge') mean_convolved = mean_convolved.astype('uint8') img_show2(img, mean_convolved)
平均化フィルタによる平滑化の結果, ブラーがかかっていることがわかります。
3.1.2 opencvによる実装
OpenCVでは, cv2.filter2D 関数を用いることで任意カーネルによる畳み込みができ, cv2.blur 関数を用いることで平均化フィルタを利用できるようです。
img = cv2.cvtColor(cv2.imread('img_src/wizdomiot.png',1), cv2.COLOR_BGR2RGB) blur = cv2.blur(img,(29,29)) img_show2(img, blur)
手実装したものと速度比較してみます。
%%timeit blur = cv2.blur(img,(29,29))
3.55 ms ± 235 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%%timeit N = 29 kernel = np.ones([N,N]) / (N*N) mean_convolved = np.zeros(img.shape) for (i,ch) in enumerate(['R','G','B']): mean_convolved[:,:,i] = convolve2d_1ch(img[:,:,i], kernel, padding='edge') mean_convolved = mean_convolved.astype('uint8')
4.07 s ± 507 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
N=29 で 1000 倍ほどの速度差が出ました。OpenCVはすごいですね。
また, N=5 などの小さい数で比較してみると,
%%timeit blur = cv2.blur(img,(5,5))
2.68 ms ± 300 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%%timeit N = 5 kernel = np.ones([N,N]) / (N*N) mean_convolved = np.zeros(img.shape) for (i,ch) in enumerate(['R','G','B']): mean_convolved[:,:,i] = convolve2d_1ch(img[:,:,i], kernel, padding='edge') mean_convolved = mean_convolved.astype('uint8')
284 ms ± 6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
このように, OpenCVでは N=5 と N=29 で速度差があまりない一方で, 畳み込み演算を素朴に実装した方は N に比例して大きくなってしまっていることがわかります。
畳み込み演算の計算を考えると, こちらの時間計算量は $O(N^2)$ といったところでしょうか。
3.2. ガウシアンフィルタ
平均化フィルタでは, 目的のピクセルとその周辺のピクセルを同じ重みづけで計算をしていたことになりますが, 目的のピクセルに近いピクセルの比重を大きくすることも考えられます。
そこで, 目的のピクセルを中心とした ガウス分布(正規分布) に従って重み付き平均を計算するフィルタのことを, ガウシアンフィルタ と呼びます。
今回使うカーネルは二次元なので, 以下のような二次元ガウス分布になります。
平均:0, 分散: $ σ^2$ の二次元正規分布:
$ h_g(x,y) = \dfrac{1}{2 \pi σ^{2}} \exp \bigg( - \dfrac{ x^{2}+y^{2}}{2 σ^{2} } \bigg) $
とりあえずこれをそのまま関数化して,
# 平均 0, 分散 sigma^2 の二次元ガウス分布 def norm2d(x,y,sigma): Z = np.exp(-(x**2 + y**2) / (2 * sigma**2)) / (2 * np.pi * sigma**2) return Z
プロットしてみると,
sns.set_style('ticks') fig = plt.figure(figsize=(10,8)) ax = fig.add_subplot(111, projection='3d') x = y = np.arange(-3,3,0.1) X,Y = np.meshgrid(x,y) ax.plot_surface(X,Y,norm2d(X,Y,1),cmap='coolwarm') plt.show()
このような感じです。二次元になっただけの普通の正規分布です。
これを平均化カーネルの重み付けに利用するのですが,
「フィルタサイズ (2w + 1)(2w + 1) の場合，σ = w/2 とするのが一つの目安」 http://www.clg.niigata-u.ac.jp/~lee/jyugyou/img_processing/medical_image_processing_03_press.pdf
というような記述があったので, ひとまずこれに従ってみることにします。
以上を踏まえて, ガウシアンフィルタのためのカーネル行列を作成する関数を以下に実装しました。
def gaussian_kernel(size): if size%2==0: print('kernel size should be odd') return sigma = (size-1)/2 # [0,size]→[-sigma, sigma] にずらす x = y = np.arange(0,size) - sigma X,Y = np.meshgrid(x,y) mat = norm2d(X,Y,sigma) # 総和が1になるように kernel = mat / np.sum(mat) return kernel
以下のようになります。
gaussian_kernel(3)
array([[ 0.07511361, 0.1238414 , 0.07511361],
[ 0.1238414 , 0.20417996, 0.1238414 ],
[ 0.07511361, 0.1238414 , 0.07511361]])
gaussian_kernel(5)
array([[ 0.02324684, 0.03382395, 0.03832756, 0.03382395, 0.02324684],
[ 0.03382395, 0.04921356, 0.05576627, 0.04921356, 0.03382395],
[ 0.03832756, 0.05576627, 0.06319146, 0.05576627, 0.03832756],
[ 0.03382395, 0.04921356, 0.05576627, 0.04921356, 0.03382395],
[ 0.02324684, 0.03382395, 0.03832756, 0.03382395, 0.02324684]])
中央が一番大きく, 周辺に向かうにつれて割合が減衰していることがわかると思います。
カーネルが用意できたので, 早速畳み込んでみましょう。
gauss_convolved = np.zeros(img.shape) for (i,ch) in enumerate(['R','G','B']): gauss_convolved[:,:,i] = convolve2d_1ch(img[:,:,i], gaussian_kernel(29), padding='edge') gauss_convolved = gauss_convolved.astype('uint8') img_show2(img, gauss_convolved)
かかりました。先ほどの平均化フィルタとの違いがよくわからないので, 並べてみます。
img_show2(mean_convolved, gauss_convolved,['mean','gaussian'])
やや微妙な差ですが, ガウシアンフィルタのほうがなめらかな結果になっています。
3.3. 応用 : 特定方向の平滑化
上記の応用として, 特定方向のみの平滑化, いわゆるブラー(方向)があります。これは, 特定方向のみの平均を取ることによって実現できます。
# 5*5 単位行列を5で割る np.eye(5)/5
array([[ 0.2, 0. , 0. , 0. , 0. ],
[ 0. , 0.2, 0. , 0. , 0. ],
[ 0. , 0. , 0.2, 0. , 0. ],
[ 0. , 0. , 0. , 0.2, 0. ],
[ 0. , 0. , 0. , 0. , 0.2]])
例えば, 上記のようなカーネルを用いることで, 斜め方向の平滑化が可能になります。
N=39 kernel = np.eye(N)/N convolved_1 = np.zeros(img.shape) for (i,ch) in enumerate(['R','G','B']): convolved_1[:,:,i] = convolve2d_1ch(img[:,:,i], kernel, padding='edge') convolved_1 = convolved_1.astype('uint8') img_show2(img, convolved_1)
斜めにかかりました。しかし, あまり仕上がりが綺麗ではありません。
そこで, 先ほどと同じように, 正規分布に従って重み付けをしてみることにします。シンプルに以下のようにします。
mat = gaussian_kernel(5)*np.eye(5) mat / np.sum(mat)
array([[ 0.11170336, 0. , 0. , 0. , 0. ],
[ 0. , 0.23647602, 0. , 0. , 0. ],
[ 0. , 0. , 0.30364122, 0. , 0. ],
[ 0. , 0. , 0. , 0.23647602, 0. ],
[ 0. , 0. , 0. , 0. , 0.11170336]])
N=39 mat = gaussian_kernel(N)*np.eye(N) kernel = mat / np.sum(mat) convolved_2 = np.zeros(img.shape) for (i,ch) in enumerate(['R','G','B']): convolved_2[:,:,i] = convolve2d_1ch(img[:,:,i], kernel, padding='edge') convolved_2 = convolved_2.astype('uint8') img_show2(img, convolved_2)
img_show2(convolved_1,convolved_2,['mean','gaussian'])
これは右のガウス分布重み付き平均のほうが綺麗だと思うのですが, いかがでしょうか。。。
ここで, 任意方向の平滑化を作ってみたいと思います。カーネルの作り方は色々ありえると思いますが, ここでは cv2.line() 関数を用いてみます。
具体的には, cv2.line() では, 始点と終点の座標を指定することで以下のように「線を引く」ことができるので,
N=7 bg = np.zeros([N,N]) bg_lined = cv2.line(bg,(2,0),(4,7),1,1) print(bg_lined)
[[ 0. 0. 1. 0. 0. 0. 0.]
[ 0. 0. 1. 0. 0. 0. 0.]
[ 0. 0. 0. 1. 0. 0. 0.]
[ 0. 0. 0. 1. 0. 0. 0.]
[ 0. 0. 0. 1. 0. 0. 0.]
[ 0. 0. 0. 0. 1. 0. 0.]
[ 0. 0. 0. 0. 1. 0. 0.]]
これを重み付き平均の係数の形にしてあげればよい,ということです。
gauss_lined = bg_lined * gaussian_kernel(N) kernel = gauss_lined / np.sum(gauss_lined) print(kernel)
[[ 0. 0. 0.10329792 0. 0. 0. 0. ]
[ 0. 0. 0.13637317 0. 0. 0. 0. ]
[ 0. 0. 0. 0.17030948 0. 0. 0. ]
[ 0. 0. 0. 0.18003887 0. 0. 0. ]
[ 0. 0. 0. 0.17030948 0. 0. 0. ]
[ 0. 0. 0. 0. 0.13637317 0. 0. ]
[ 0. 0. 0. 0. 0.10329792 0. 0. ]]
これを, 任意の角度の入力に対して行えるようにしたカーネル作成関数を実装します。
始点と終点の座標は行列の外側でもかまわないことに注意すると, 以下のようになるかと思います。
def directional_kernel(size,angle): if size%2==0: print('kernel size should be odd') return angle = angle % 180 radian = (angle/180) * np.pi x1 = int(np.cos(radian)*(size/2) + (size/2)) y1 = int(np.sin(radian)*(size/2) + (size/2)) x2 = int(np.cos(radian + np.pi)*(size/2) + (size/2)) y2 = int(np.sin(radian + np.pi)*(size/2) + (size/2)) bg = np.zeros([size,size]) bg_lined = cv2.line(bg,(x1,y1),(x2,y2),1,1) gauss_lined = bg_lined * gaussian_kernel(size) kernel = gauss_lined / np.sum(gauss_lined) return kernel
実行してみます。
directional_kernel(5,0)
array([[ 0. , 0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. ],
[ 0.15246914, 0.2218413 , 0.25137912, 0.2218413 , 0.15246914],
[ 0. , 0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. ]])
directional_kernel(5,30)
array([[ 0. , 0. , 0. , 0. , 0. ],
[ 0.14343172, 0.20869192, 0. , 0. , 0. ],
[ 0. , 0. , 0.26796573, 0.23647892, 0. ],
[ 0. , 0. , 0. , 0. , 0.14343172],
[ 0. , 0. , 0. , 0. , 0. ]])
directional_kernel(5,60)
array([[ 0. , 0.14343172, 0. , 0. , 0. ],
[ 0. , 0.20869192, 0. , 0. , 0. ],
[ 0. , 0. , 0.26796573, 0. , 0. ],
[ 0. , 0. , 0.23647892, 0. , 0. ],
[ 0. , 0. , 0. , 0.14343172, 0. ]])
このような形になりました。実際に使ってみると,
def directinal_blur(img, N, angle): convolved = np.zeros(img.shape) for (i,ch) in enumerate(['R','G','B']): convolved[:,:,i] = convolve2d_1ch(img[:,:,i], directional_kernel(N,angle), padding='edge') convolved = convolved.astype('uint8') img_show2(img, convolved)
directinal_blur(img,N=39,angle=0)
directinal_blur(img,N=39,angle=30)
directinal_blur(img,N=39,angle=60)
directinal_blur(img,N=39,angle=90)
0,30,60,90度の方向にブラーがかかっていることがわかります。
4. おわりに
今回は, 畳み込みといくつかの平滑化フィルタを実装することで, 空間フィルタリングに関する初歩的な理解を得ることができました。
次回は, メディアンなどのエッジを保存する平滑化フィルタについて扱おうかと思います。
