ガウシアンフィルタ

ガウシアンフィルタとはなにか

画像処理における ガウシアンフィルタ(Gaussian filter) とは、 平滑化フィルタ(smoothing filter) の一つであり、画像に含まれるノイズを軽減するのに使われる。

ガウス分布を近似した カーネル(kernel) を利用して畳み込み演算をすることで、フィルタリングを行う。

2次元ガウス分布

画像のガウシアンフィルタでは、2次元ガウス分布を利用するが、特別考慮することがなければ、x方向とy方向は独立で、同じ分散に従っているようなガウス分布を考える。

つまり、xとyに相関があったり分散が異なるような一般的なガウス分布を考える必要はない。

まずそういった2次元ガウス分布の確率密度関数を表示してみよう。

x方向の周辺分布が、標準偏差σの1次元ガウス分布

$$ g_x(x) = \frac{1}{\sqrt{2 \pi \sigma ^ 2}} \exp{(-\frac{x^2}{2 \sigma^2})} $$

y方向の周辺分布も、標準偏差σの1次元ガウス分布

$$ g_y(y) = \frac{1}{\sqrt{2 \pi \sigma ^ 2}} \exp{(-\frac{y^2}{2 \sigma^2})} $$

である。 これらが独立なので、求める確率密度関数は、

$$ G(x, y) = g_x(x) g_y(y) = \frac{1}{2 \pi \sigma ^ 2} \exp{(-\frac{x^2 + y^2}{2 \sigma^2})} $$

である。

2次元ガウス分布の近似

この分布$G(x, y)$は無限に広がっている。

無限に連続的に広がっているものを有限で離散的なピクセルの集合である画像と畳み込み計算するのを考えるのは難しいので、 この$G(x, y)$を奇数かける奇数の有限の離散的な格子に近似したものがカーネル(kernel)である。

カーネルの具体例としては、よく次のようなものが使われる。

これらのカーネルはどのように$G(x, y)$を近似したのだろうか？

実は、各格子に対して$G(x, y)$を計算して、全体の和が1になるように正規化したのである。

たとえば、左の3x3の格子の場合、$σ=0.85$とすると

$$ G(0, 0) = \frac{1}{2 \pi \sigma ^ 2} \exp{(-\frac{0}{2 \sigma^2})} = \frac{1}{2 \pi \sigma ^ 2} \fallingdotseq 0.22 $$

$$ G(1, 0) = G(-1, 0) = G(0, 1) = G(0, -1) = \frac{1}{2 \pi \sigma ^ 2} \exp{(-\frac{1}{2 \sigma^2})} \fallingdotseq 0.11 $$

$$ G(1, 1) = G(-1, 1) = G(1, -1) = G(-1, -1) = \frac{1}{2 \pi \sigma ^ 2} \exp{(-\frac{2}{2 \sigma^2})} \fallingdotseq 0.055 $$

となり、これらの和は

$$ S = \sum_{-1 \le x \le 1, -1 \le y \le 1}{G(x, y)} = 0.22 + 0.11 \cdot 4 + 0.055 \cdot 4 = 0.88 $$

となる。そこでカーネルの各格子の値は

$$ K(0, 0) = \frac{G(0, 0)}{S} = \frac{0.22}{0.88} = \frac{4}{16} $$

$$ K(1, 0) = K(-1, 0) = K(0, 1) = K(0, -1) = \frac{G(1, 0)}{S} = \frac{0.11}{0.88} = \frac{2}{16} $$

$$ K(1, 1) = K(-1, 1) = K(1, -1) = K(-1, -1) = \frac{G(1, 1)}{S} = \frac{0.055}{0.88} = \frac{1}{16} $$

と計算される。

各$G(x, y)$を計算をするとき、最後に結局正規化をするので、$\frac{1}{2 \pi \sigma^2}$は乗じても乗じなくても良い。

乗じない場合を$G’(x, y)$とおけば、

$$ G’(x, y) = \exp{(-\frac{x^2 + y^2}{2 \sigma^2})} $$

また、

$$ r = \exp{(-\frac{1}{2 \sigma^2})} $$

とおくと、

$$ G’(x, y) = r ^ {x^2 + y^2} $$

と非常に簡潔になる。

つまり、カーネルを計算するうえでは、まず$\sigma$から$r$を求め、カーネルの中心からの距離の2乗分だけ$r$を累乗した格子を書く。

これを正規化すると、カーネルが得られる。

逆にカーネルが与えられたとき、そのベースとなるガウス分布の$\sigma$を求めるには、中心と中心の一つ横の格子の値の比をとって、$r$とし、

$$ \sigma = \sqrt{-\frac{1}{2 \log{r}}} $$

によって求めれば良い。

たとえば、上に示した5x5の格子の場合、

$$ r = \frac{\frac{24}{256}}{\frac{36}{256}} = \frac{24}{36} = \frac{2}{3} $$

であるから、

$$ \sigma = \sqrt{-\frac{1}{2 \log{r}}} = \sqrt{-\frac{1}{2 \log{\frac{2}{3}}}} \fallingdotseq 1.11 $$

となる。

以上の議論を踏まえて、numpyを使用してガウシアンフィルタのカーネルを計算するコード例は次のようなものである。

import numpy as np


def make_gaussian_kernel(width: int, height: int, sigma: float) -> np.ndarray
    # kernelの中心と1個横のセルの値の比
    r = np.exp(-1 / (2 * sigma**2))

    # kernelの中心index
    mw = width // 2
    mh = height // 2

    # kernelの中心を原点とした座標が入った配列
    xs = np.arange(-mw, mw + 1)
    ys = np.arange(-mh, mh + 1)
    x, y = np.meshgrid(xs, ys)

    # rを絶対値の2乗分累乗すると、ベースとなる配列が求まる
    g = r**(x**2 + y**2)

    # 正規化
    return g / np.sum(g)

畳み込み

こうしてできたカーネルを、元の画像に畳み込むことでガウシアンフィルタとなる。

畳み込み演算(convolution) は、「カーネルに従って周りの画素値を取り込むこと」、と考えればよい。

元の画像の画素値を$I(x, y)$として、カーネルを$K(x, y)$とすれば、出力画像の画素値$J(x, y)$は次のようになる。

$$ J(x, y) = K * I = \sum_{u=-w}^{w} \sum_{v=-h}^{h} K(u, v) I(x + u, y + v) $$

ここで、$w$はカーネルの幅の半分、$h$はカーネルの高さの半分である。

各点の周りで、カーネルとの内積を取る (要素ごとにかけて足す) ことをしているわけである。

畳み込みを行う際は、元の画像の端っこの処理についていくつかのバリエーションが考えられる。

たとえば、scipyのconvolve2dには、いくつかのモードが用意されている。

https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.convolve2d.html

フィルタ例

以下は、実際に3x3、$\sigma = 0.85$のカーネルを使用して、ガウシアンフィルタを行った例である。

元の画像が平滑化されていることがわかる。

エッジ部分もぼやけてしまうため、ノイズの性質によって適切なフィルタを選択する必要はある。

元画像

フィルタ後画像