# 画像処理 エッジ&コーナー検出のための空間フィルタリング(1)

記事名に(1)をつけたのは、まだ実装できていないノンローカルミーンフィルタを後ほど実装したいと考えているためです。

・追記(2022/02/13)
　ノンローカルミーンフィルタについての記事を書きました。

当ページでは、理解を深めるために、OpenCVに頼らないフィルタの実装方法を掲載しています。
コードに間違い等があればご指摘していただけると幸いです。

空間フィルタリングとは

空間フィルタリングとは、画像処理を使ってエッジやコーナー検出を行う前に画像にあらかじめ施しておく処理のことです。
これを行う理由には、写真を撮った時に発生するノイズによる影響をできるだけ防ぐために行われます。

夜景を撮る時など、低光量であるときにノイズの影響はより大きく鮮明に現れるようになります。
夜景を撮ったら以下の写真のように暗い部分が鮮明な暗さではなくなったという経験はあるのではないでしょうか?

これはカメラの不良などによるノイズが含まれてしまった結果です。
(ちょっとわかりにくいですが)

このノイズは性能が悪いカメラほど顕著に現れます。

だったら、性能が良いカメラで撮ったものを使えばいいじゃないか、となるかもしれませんが、実際にはドライブレコーダーや防犯カメラなど、画質は劣るが安価で量産に向いており大雑把に撮影することで十分とされるカメラを対象とした場合は、フィルタリングが不要なくらいに画質を求めるのは難しいと考えられます。
ゆえに、汎用性を高めるためには、空間フィルタリングが重要となるのではないかと個人的には考えています。

空間フィルタリングの適用例

空間フィルタリングは、その注目している画素値とその近傍の画素値を入力値として、出力する画素に対応する画素値を計算する方法です。周りの値を用いて何かしらの計算を行うことで、ノイズの影響を軽減することができます。

例えば、以下のような白黒画像を検討します。
マスは一つの画素であり、それらは0から255までの値で明るさを表しており、大きければ大きいほど明るいと考えてください。

試しに、注目する画素1つとその周りの画素8つの値の平均を取るようなフィルタを検討します。その場合だと、各画素は以下のように計算することができます。


そしてフィルタをかけた後の画像は以下のようになります。

すると、黒しかない画像のなかに白いノイズ点が含まれていた場合でも、空間フィルタリングをかけることによって、この白い点をノイズとして平滑化を行い、周りと同じような色にして目立たなくしてくれることがわかります。
これが主なフィルタリングの役割です。

ただ、このフィルタリング効果が強くなってしまうと、色の濃淡がはっきりとすべき箇所にまで影響を及ぼしてしまう(ピンボケが酷い画像のようになる)ため、適切な範囲で施す必要があります。

空間フィルタリングの数学的定義

空間フィルタリングを数学的定義に基づいて書くと以下のようになります。

入力画素の値と出力画素の値をそれぞれ $I,J$ とします。
座標 $(x, y)$ 、近傍領域の範囲を $(2N+1)\times(2N+1)$ として、空間フィルタリングをかけた時の画素の計算方法は以下のように表示できます。(ただし、$N$ は1以上の整数)

J(x, y) = \sum_{j=-N}^{N} \sum_{i=-N}^{N} F(i, j)I(x+i, y+j)

ここでの $F$ はカーネルと呼ばれ、注目する画素に対してどのような計算方法を適用するかを行列で定義します。$N$ の値はカーネルの行数と列数によって決まります。

また、カーネル $F$ を計算する方法として、相関と畳み込みの2種類があります。
相関とは、カーネルと入力画素の行と列が一致しているもの同士を計算する方法で、畳み込みは行と列を逆に並び替えたカーネルと入力画素を計算する方法です。文字で説明してもわかりにくいので、以下の画像を元に理解されると良いと思います。

数式で表すと以下の通りです。
画素 $I$ とカーネル $F$ がそれぞれ

I = 
\begin{pmatrix}
a & b & c\\
d & e & f\\
g & h & i
\end{pmatrix}
,
\hspace{10px}

F = 
\begin{pmatrix}
r & s & t\\
u & v & w\\
x & y & z
\end{pmatrix}

である時の計算方法は以下の通りです。

相関	畳み込み
$J = ar + bs + ct + du + ev + fw + gx + hy + iz$	$J = az + by + cx + dw + ev + fu + gt + hs + ai$

数式で書くと複雑になってしまいますが、先程の実際に計算してみた例では、$N = 1$としてそれぞれの画素の計算を行うことと同じ意味となります。

単純な計算でも、数式だと難しいことをやっているように表記できますが、実際のところは大したことをやっていないというのが現実です。

空間フィルタリングの種類

今回はさまざまな空間フィルタリングをPythonで実装していこうと思います。
なお、適用する画像は以下のランドクルーザーの写真を用います。
また、結果をわかりやすくするために、あえてノイズを付加しております。

(画像はこちらの「Response」さんのランドクルーザーの記事から引用させていただきました)

1.平均化フィルタ, 中央値フィルタ

平均化フィルタは、先ほど例に示したような、決めた範囲について平均値を画素値に適用するフィルタののことであり、中央値フィルタは、決めた範囲の中央値を画素値に適用するフィルタのことです。

これらは、先程の例のように単純な計算をすれば実装できますので、Pythonでの実装は割愛します。

2.ガウシアンフィルタ

ガウシアンフィルタでは、注目する画素から近いと大きい重み付けを、遠いと小さい重み付けをするようなカーネルを用いたフィルタリングを行います。計算式は以下のように表され、分散$\sigma$の値を変えることで、注目画素とその周りの画素との重みの差異を変えることが可能です。式は正規分布に似ています。

F(x, y) = \frac{1}{2 \pi \sigma ^ 2} \exp \Bigl( -\frac{x^2 + y^2}{2\sigma ^ 2}\Bigr)

この式の$x$と$y$に対する$F(x, y)$の分布は以下のグラフのようになります。中央であるほど値が大きく、離れるにつれて小さくなっていることが確認できます。

この分布をカーネルに適用してフィルタリングを行います。

Pythonでの実装

import cv2
import numpy as np

# 必要な関数の定義
# カーネル関数
def kernel_function(x, y, sigma):
    z = (1 / (2 * np.pi * sigma ** 2)) * np.exp(-(x ** 2 + y ** 2) / (2 * sigma ** 2))
    return z

# カーネルを作成する関数
def create_kernel(size=(3, 3), sigma=1):
    kernel = np.zeros(size)
    kerner_x, kernel_y = int((size[1] - 1) / 2), int((size[0] - 1) / 2)
    for x in range(size[1]):
        for y in range(size[0]):
            kernel[y][x] = kernel_function(x - kerner_x, y - kernel_y, sigma)
    # カーネルの規格化
    kernel_norm = kernel / np.sum(kernel)
    return kernel_norm

# ガウシアンフィルタ
def gaussian_filter(img, kernel):
    # カーネルサイズを取得
    line, column = kernel.shape
    # フィルタをかける画像の高さと幅を取得
    height, width = img.shape

    # 畳み込み演算をしない領域の幅を指定
    height_ignore = int((line - 1) / 2)
    width_ignore = int((column - 1) / 2)

    # 出力画像用の配列
    img_filter = img.copy()

    # フィルタリングの計算
    for y in range(height_ignore, height - height_ignore):
        for x in range(width_ignore, width - width_ignore):
            img_filter[y][x] = np.sum(img[y - height_ignore: y + height_ignore + 1, x - width_ignore: x + width_ignore + 1] * kernel)
    return img_filter

上記では必要な関数を定義しています。
kernel_functionでは、先程の$F(x, y)$の式を計算する関数です。
create_kernelでは、カーネルの大きさと分散の値に従ってカーネルを計算します。
gaussian_filterでは、上記2つの関数を使って実際に画像にフィルタをかける関数です。

# 入力画像を白黒画像として読み込み
img_gray = cv2.imread("rankuru_noise.jpg", cv2.IMREAD_GRAYSCALE)

# カーネルを定義
# sizeにはカーネルを適用する計算範囲を、sigmaには分散を渡す
kernel = create_kernel(size=(5, 5), sigma=1)

# ガウシアンフィルタを適用
img_filter = gaussian_filter(img_gray, kernel)

先程の3つの関数を使ってフィルタをかけると、フィルタを適用した画像のデータが得られます。
フィルタ適用前と適用後にわけて比較を行いましょう。可視化にはmatplotlibライブラリのimshow関数を用います。

# 元の画像とフィルタ適用後の画像を比較
import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(20, 30))
ax1.imshow(img_gray, cmap='gray')
ax2.imshow(img_filter, cmap='gray')
ax1.set_xticks([]); ax1.set_yticks([])
ax2.set_xticks([]); ax2.set_yticks([])
plt.show()

結果は以下のようになると思います。

微妙にフィルタ効果があるんですが、わかりにくいですね...
使用する画像が悪かったかもですが...
あとよーく見ると、フィルタをかけた後に車の輪郭がぼやけていることがわかると思います。
ガウシアンフィルタでは、平均化フィルタよりぼやけ具合は抑えられますが、それでもまだ少しはぼやけてしまうという問題点があります。

このぼやけを無くす、つまり、必要なところだけ平滑化を行い、輪郭をしっかり残そうとして考案された方法が、次に紹介するバイラテラルフィルタです。

3.バイラテラルフィルタ

バイラテラルフィルタとは、ガウシアンフィルタの考え方に加えて、画素値に依存する重み付けを行うフィルタのことを言います。
注目画素との値の差異が大きければ重みを小さくし、逆に、差異が小さければ、重みを大きくすることによって、ガウシアンフィルタでは問題となっていた輪郭ボケを抑えることができるようになります。

計算式は以下の通りになっています。$I(x,y)$はフィルタをかける前の画素値を、$\tilde{I}(x,y)$はフィルタをかけた後の画素値を、$W(x,y,i,j)$は適用するカーネルを表しています。$\exp$の指数部分の第二項が新たに追加された部分になります。

\tilde{I}(x,y) = \frac{\sum_{i=-N}^{N} \sum_{j=-N}^{N} W(x,y,i,j)I(x,y)}{\sum_{i=-N}^{N} \sum_{j=-N}^{N} W(x,y,i,j)}

W(x,y,i,j) = \exp \Biggl(-\frac{i^2+j^2}{2\sigma_d^2} - \frac{||I(x,y) - I(x+i, y+j)||^2}{2\sigma_r^2}\Biggr)

特徴としては、一つの画素の計算毎にカーネルが毎度変わることです。画素値に基づいてカーネルを計算するため、毎度カーネルを計算し直す必要があります。

Pythonでの実装

import cv2
import numpy as np

# カーネル関数
def kernel_function(x, y, I, I_dash, sigma_d, sigma_r):
    z = np.exp(-((x ** 2 + y ** 2) / (2 * sigma_d ** 2)) - (((int(I) - int(I_dash)) ** 2) / (2 * sigma_r ** 2)))
    return z

# カーネルを作成する関数
def create_kernel(matrix, sigma_d, sigma_r):
    size = matrix.shape
    kernel = np.zeros(size)
    kernel_x, kernel_y = int((size[1] - 1) / 2), int((size[0] - 1) / 2)
    kernel_center = matrix[kernel_y][kernel_x]
    for x in range(size[1]):
        for y in range(size[0]):
            kernel[y][x] = kernel_function(x - kernel_x, y - kernel_y, kernel_center, matrix[y][x], sigma_d, sigma_r)
    return kernel

# バイラテラルフィルタの演算を行う関数
def bilateral_function(matrix, kernel):
    I_bar = np.sum(matrix * kernel) / np.sum(kernel)
    return I_bar

# バイラテラルフィルタ
def bilateral_filter(img, size=(3, 3), sigma_d=1, sigma_r=1):
    # フィルタをかける画像の高さと幅を取得
    height, width = img.shape

    # 畳み込み演算をしない領域の幅を指定
    height_ignore, width_ignore = int((size[0] - 1) / 2), int((size[1] - 1) / 2)

    # 出力画像用の配列
    img_filter = img.copy()

    # フィルタリングの計算
    for y in range(height_ignore, height - height_ignore):
        for x in range(width_ignore, width - width_ignore):
            # カーネルを作成
            matrix = img[y - height_ignore: y + height_ignore + 1, x - width_ignore: x + width_ignore + 1]
            kernel = create_kernel(matrix, sigma_d, sigma_r)
            # カーネルを用いて計算
            img_filter[y][x] = bilateral_function(matrix, kernel)
    return img_filter

ガウシアンフィルタより少し複雑になっていますが、大まかな流れは同じです。
kernel_functionでは、先程の$W(x, y, i, j)$の式を計算する関数です。
create_kernelとbilateral_functionの2つの関数を使って、カーネルの大きさと分散の値に従ってカーネルを計算します。

# 入力画像を読み込み
img_gray = cv2.imread("rankuru_noise.jpg", cv2.IMREAD_GRAYSCALE)

# バイラテラルフィルタを適用
# sizeにはカーネルを適用する計算範囲を、sigma_dとsigma_rにはそれぞれ距離に対する分散と画素値の差異に対する分散を渡す
img_filter = bilateral_filter(img_gray, size=(15, 15), sigma_d=20, sigma_r=20)

先程の関数を使ってフィルタをかけると、フィルタを適用した画像のデータが得られます。
matplotlibライブラリを使って可視化してみましょう。

import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(20, 30))
ax1.imshow(img_gray, cmap='gray')
ax2.imshow(img_filter, cmap='gray')
ax1.set_xticks([]); ax1.set_yticks([])
ax2.set_xticks([]); ax2.set_yticks([])
plt.show()

結果から、フィルタによる平滑化ができており、かつ、輪郭がぼやけずにしっかりと残っていることがわかります。

ただ、計算効率が非常に悪いためか、フィルタ適用に5分もかかってしまいました。
計算を早めるためには工夫や改修が必要になるかと思われますが、そもそもOpenCVライブラリを使えば速度は解決するため、ここでは計算ができればよしということで。

まとめ

空間フィルタリングについてと、その代表例である「ガウシアンフィルタ」と「バイラテラルフィルタ」のPythonによる自力実装を行いました。
「ノンローカルミーンフィルタ」などもありますが、これについては実装が出来次第記事にしようと思います。
((2022/02/13) 冒頭にある通り、ノンローカルミーンフィルタについての記事は投稿済みです)

以上になります。