ガボールフィルタの概要

信号処理の技法には、多数のフィルタが存在するが、ここで科学的な応用のきわめて多いガボールフィルタ(Gabor filter)に着目する。我々の目的は、画像から、多種多様な情報を「特徴量」として取り出し、最終的に控えている予測課題の「入力」として使うことである。詳しい紹介は専門書に委ねるとして、ここでは要点のみ述べることにする。

一言でいえば、ガボールフィルタとは、ガウス関数で変調される正弦波である。英語の呼称を使って定式化すると、

\begin{align} (\text{Filter output}) = (\text{Sinusoid carrier}) \times (\text{Gaussian envelope}) \end{align}

という形式になる。その「次元数」によって1D、2D、3Dと様々なバージョンが文献の上では存在する。その「カタチ」をイメージする上で、たとえばJones et al. (1987)では明快な図が多数ある。応用事例をいくつかピックアップしてみると、たとえば、動物の受容野のモデリング(DeAngelis et al., 1993)、テキスチャー・質感の判定(scikit-image, "Gabor filter banks for texture classification")、顔認識(九工大のBIS研より)、など、実に多様な応用例がある。我々の実務として重要になってくるのは、下記の事項である:

  • 周波数のパラメータを調整することで、フィルタが「効く」解像度が変わる。

  • ガウス関数の半径を大きく・小さくすることで、画像内で効く範囲が拡大・縮小される。

  • もっとも単純なガボールフィルタは1Dで、時間の関数として信号を捉えていることから「temporal」と呼ぶ。2次元以上だと、空間「spatial」や時空「spatio-temporal」のバリエーションがある。

フィルタを信号に適用する方法はいくつもあるが、もっとも標準的なのは、畳み込み(convolution)という操作である。対象信号とフィルタ(これもまた信号である)の定義域が共通しているので、畳み込んで新しい信号を造ることができる。表記についてだが、任意の信号$f$と$g$をとって、これらの畳み込みを$f \ast g$と記す。1次元の場合は下記のように定義される。

\begin{align} (f \ast g)(t) = \int_{\mathbb{R}} f(t)g(t-\tau) \, d\tau \end{align}

2次元もまったく同様に定義される。

\begin{align} (f \ast g)(x,y) = \int_{\mathbb{R}^{2}} f(x,y)g(x-\alpha,y-\beta) \, d\alpha\,d\beta. \end{align}

呼称はいくらでもあり得るが、もし$f$が関心のある入力(ここでは画像内の画素値を返す)で、$g$が我々が設計したフィルタであるとすると、フィルタを応用した結果として得られる応答が$f \ast g$である。

1次元のガボールフィルタ

概要と良いとして、これから自分でフィルタを設計・実装しようと思えば、もう少し具体的な情報が必要だ。まず、基本的な表記とフィルタの働きを確認するために、1Dの「temporal」ガボールフィルタを作ってみよう。書き方として、

\begin{align} G(t) = f(t) \, s(t) \end{align}

とするが、これらのパーツは下記のように定義される。

\begin{align} s(t) & = \exp \left(i (2\pi ut + \phi) \right)\\ & = \cos(2\pi ut + \phi) + i\sin(2\pi ut + \phi)\\ f(t) & = A \exp\left( -\frac{t^2}{\sigma^2} \right). \end{align}

正弦波とガウス関数の効き目をつかさどるパラメータについて:

  • $u$: 周波数。つまり、正弦波が時間1単位あたり、何周するか示す値。例として、$u=1$なら、周期が1(時間、分、秒など)である。周波数を$u=2$と倍にすると、周期が半減して1/2になる(時間1単位あたり、2周もするから)。反対に$u=1/2$ならば、一周するには時間が2単位必要になる。以下同様。
  • $A$: 振幅。フィルタの取る最大値を定める。
  • $\phi$: 位相。定義域において、どこで(時間ならば、いつ)大きく・小さくなるか決める。
  • $\sigma$: ガウス関数の範囲を定める。 $\sigma > 0$が小さいほど、効き目を持つ時間が短くなる。

数式はもう良いとして、そろそろ自分で実装して、ガボールフィルタを動かしてみよう。

In [1]:
import math
import numpy as np


def G_carrier_real(t, freq, phase):
    '''
    Real part of the carrier.
    '''
    topass = 2 * math.pi * freq * t + phase
    out = np.cos(topass)
    return out

def G_carrier_imag(t, freq, phase):
    '''
    Imaginary part of the carrier.
    '''
    topass = 2 * math.pi * freq * t + phase
    out = np.sin(topass)
    return out


def G_envelope(t, amp, sdev):
    '''
    The impact of the filter is controlled by a Gaussian function.
    '''
    out = amp * np.exp( (-(t/sdev)**2) )
    return out

フィルタ自体は、複数のパラメータによって決められる。下記のdictでそれらの値を格納する。

In [2]:
myparas = {"amp": 1/(2*math.pi),
           "sdev": 1,
           "freq": 1,
           "phase": 0}

パラメータを引数として、カスタマイズされたフィルタを下記のように造る。

In [3]:
def G_fil_real(t, paras):
    '''
    Custom-built filter response (real part).
    Assumes that t is an array of temporal inputs.
    '''
    carrier = G_carrier_real(t=t, freq=paras["freq"], phase=paras["phase"])
    envelope = G_envelope(t=t, amp=paras["amp"], sdev=paras["sdev"])
    out = carrier * envelope
    return out

def G_fil_imag(t, paras):
    '''
    Custom-built filter response (imaginary part).
    Assumes that t is an array of temporal inputs.
    '''
    carrier = G_carrier_imag(t=t, freq=paras["freq"], phase=paras["phase"])
    envelope = G_envelope(t=t, amp=paras["amp"], sdev=paras["sdev"])
    out = carrier * envelope
    return out

上記から明らかなように、「正弦波×ガウス関数」という形式が守られている。各要素を可視化してみる。

In [5]:
import matplotlib.pyplot as plt

myfig = plt.figure(figsize=(16,4))
t_inputs = np.linspace(-3,3,500)

ax_carrier = myfig.add_subplot(1,3,1)
plt.title("Sinusoidal carrier")
plt.xlabel("Time (s)")
ax_envelope = myfig.add_subplot(1,3,2)
plt.title("Gaussian envelope")
plt.xlabel("Time (s)")
ax_filter = myfig.add_subplot(1,3,3)
plt.title("Gabor filter (1D)")
plt.xlabel("Time (s)")

ax_carrier.plot(t_inputs, G_carrier_real(t=t_inputs, freq=myparas["freq"], phase=myparas["phase"]))
ax_carrier.plot(t_inputs, G_carrier_imag(t=t_inputs, freq=myparas["freq"], phase=myparas["phase"]))

ax_envelope.plot(t_inputs, G_envelope(t=t_inputs, amp=myparas["amp"], sdev=myparas["sdev"]))

ax_filter.plot(t_inputs, G_fil_real(t=t_inputs, paras=myparas))
ax_filter.plot(t_inputs, G_fil_imag(t=t_inputs, paras=myparas))

plt.show()

練習問題 (A):

  1. freqという周波数パラメータに着目する。上左図で正弦波が実際に何周をしているか目視で確かめること。時間1単位あたりの回数が期待通りか。このfreqを色々な値(例:1、4、0.5、0.25など)に変えて、どう変わるか確認すること。

  2. 「標準偏差」のsdevパラメータをたとえば0.25や2に変えてみてください。どのようにフィルタの出力が変わるか。もしそのままでは見づらいのであれば、適宜にt_inputsの範囲を拡縮してみること。

  3. 位相のphaseに着目し、実部と虚部が一致するように、虚部のほうだけパラメータの値を変えてみてください($\pi/2$のズレがある)。


1Dフィルタの応用

フィルタ自体は見てきたが、それを応用することとはどういうことか、実践的に調べてみる。まず、わかりやすい例として擬似的な信号を造る。

In [6]:
para_HIGHFREQ = 5
para_LOWFREQ = 0.25

def my_signal(t):
    
    highfreq = 0.25 * np.sin((2*math.pi*para_HIGHFREQ*t))
    lowfreq = 2 * np.sin((2*math.pi*para_LOWFREQ*t))
    
    cond = (np.abs(t) <= 5)
    signal = highfreq + lowfreq
    out = np.select([cond], [signal])
    
    return out


myfig = plt.figure(figsize=(10,10))
t_inputs = np.linspace(-10, 10, 500)

plt.plot(t_inputs, my_signal(t=t_inputs))
plt.title("Some artificial signal")
plt.xlabel("Time (s)")
plt.show()

先述のとおり、フィルタを応用する場合は、対象信号との畳み込みを求めることが多い。scipisignalモジュールから、convolveという関数がここで最初の出番を迎える。

In [7]:
from scipy import signal

t_inputs = np.linspace(-10, 10, 1000)
sig_values = my_signal(t=t_inputs)

myparas = {"amp": 1/(2*math.pi),
           "sdev": 1,
           "freq": 0.25,
           "phase": 0}
fil_values_real = G_fil_real(t=t_inputs, paras=myparas)
fil_values_imag = G_fil_imag(t=t_inputs, paras=myparas)
fil_response_real = signal.convolve(sig_values, fil_values_real, mode="same")
fil_response_imag = signal.convolve(sig_values, fil_values_imag, mode="same")

myfig = plt.figure(figsize=(18,4))

ax_signal = myfig.add_subplot(1,3,1)
plt.title("Artificial signal")
plt.xlabel("Time (s)")
ax_filter = myfig.add_subplot(1,3,2)
plt.title("Gabor filter")
plt.xlabel("Time (s)")
ax_response = myfig.add_subplot(1,3,3)
plt.title("Filter response")
plt.xlabel("Time (s)")

ax_signal.plot(t_inputs, sig_values)

ax_filter.plot(t_inputs, fil_values_real)
ax_filter.plot(t_inputs, fil_values_imag)

ax_response.plot(t_inputs, fil_response_real)
ax_response.plot(t_inputs, fil_response_imag)

plt.show()

練習問題 (B):

  1. 上のコードを使って、myparasに入るfreqパラメータとして、複数の値を試してみてください。特に、para_LOWFREQpara_HIGHFREQの範囲内。値が小さいときに低周波数の成分が取り出せること、値が大きいときに高周波数の成分が取り出せることなど、実証的に確かめること。これを示すグラフをいくつか用意して保存すること。

  2. 位相と標準偏差のパラメータも、色々といじってみてください。それぞれの大小と、フィルタ出力の変化との関係について説明すること。

  3. 擬似信号であるmy_signalを、「低・高」から「低・中・高」に改造すること(つまり、2つではなく、3つの正弦波の線形和にする)。改めてフィルタをかけて、freqを適当に設定することで、それぞれの成分が抽出できるか。できるならば、その結果を図示すること。

以上、1次元のガボールフィルタをやや丁寧に見てきた。最終的に使うのは2次元のほうだが、基本的な働きはまったく同じであることから、1次元を先に見て、可視化しておくと2Dバージョンが非常に理解しやすくなるといえる。最大のポイントは、信号の「特徴」を、フィルタの設計の如何によって、好きなように取り出せることである。


2次元のガボールフィルタ

ガボールフィルタの基本的な発想、定式化、そして実務上の働きもある程度は見てきた。これを土台にして、より一般生の高いバージョンを見ていく。基本的な形式は1次元のときとまったく同様で、

\begin{align} (\text{Filter output}) = (\text{Sinusoid carrier}) \times (\text{Gaussian envelope}) \end{align}

となっている。ただし、「2次元」と呼ぶのは、定義域が「線」から「面」へと拡がるからである。実数値$t \in \mathbb{R}$ではなく、 実数ベクトル$(x,y) \in \mathbb{R}^{2}$をフィルタの引数とする。もう少し改まった形にすると、

\begin{align} G(x,y) = f(x,y) \, s(x,y) \end{align}

と書く。正弦波とガウス関数は至って自然な形で拡張される:

\begin{align} s(x,y) & = \exp \left(i (2\pi (ux + vy) + \phi) \right)\\ & = \cos(2\pi (ux + vy) + \phi) + i\sin(2\pi (ux + vy) + \phi)\\ f(x,y) & = A \exp\left( -\frac{x^2 + y^2}{\sigma^2} \right). \end{align}

今回は周波数パラメータ(角周波数)が$u$と$v$と2つあるが、それを極座標に置き換えて調整することが多い。変換は下記のとおりである。

\begin{align} u & = \omega \cos(\theta) \\ v & = \omega \sin(\theta) \end{align}

新しく出てきたパラメータが$\omega$(空間周波数)と$\theta$(フィルタの向き)である。以下に補足説明を記す。

  • $u$と$v$: 周波数。捉え方としては、画像における「ヨコ」と「タテ」方向にそれぞれは対応する。離散的なデジタル画像だと、空間単位(ピクセル)あたり請願はが何周するか示す数値である。
  • $A$: 振幅。1Dと同じ。
  • $\phi$: 位相。1Dと同じ。
  • $\sigma$: 標準偏差で、1Dと同様。ただし、この定式化だとガウス関数が円形になる。さらに一般的なやり方として、座標軸ごとに(たとえば$\sigma_{X}$と$\sigma_{Y}$を使って)範囲を定め、ガウス関数を楕円形にすることもできる。

数式はもう十分見てきたので、いよいよ2次元のガボールフィルタを作ってみることにしよう。

In [8]:
import numpy as np
import math
import matplotlib.pyplot as plt


def G2_carrier_real(x, y, freqx, freqy, phase):
    '''
    Real part of the 2-D Gabor carrier.
    '''
    topass = 2 * math.pi * (freqx*x + freqy*y) + phase
    out = np.cos(topass)
    return out


def G2_carrier_imag(x, y, freqx, freqy, phase):
    '''
    Imaginary part of the 2-D Gabor carrier.
    '''
    topass = 2 * math.pi * (freqx*x + freqy*y) + phase
    out = np.sin(topass)
    return out


def G2_envelope(x, y, amp, sdev):
    '''
    Gaussian envelope for a 2-D Gabor filter.
    We assume that it is circular (same rate of decrease in x/y directions).
    '''
    out = amp * np.exp(-(x**2+y**2)/(sdev**2))
    return out

1Dの場合と同様、いくつかのパラメータをdictにまとめ、フィルタの働きを制御する。

In [9]:
PIX_W = 128 # image width, in pixels
PIX_H = 128 # image height, in pixels
myparas = {"freqs": 4/max(PIX_W,PIX_H), # cycles per pixel
           "dir": math.pi/2, # orientation
           "amp": 1,
           "sdev": max(PIX_W,PIX_H)/5,
           "phase": 0}

ここで新たに見えるfreqsとは空間周波数(spatial frequency)で、$\omega$に対応する。パラメータをまとめて、カスタマイズされたフィルタをここでも作ってみよう。

In [10]:
def G2_fil_real(x, y, paras):
    '''
    Custom-built filter response (real part).
    '''
    # Spatial frequency in polar coordinates.
    u = paras["freqs"] * math.cos(paras["dir"])
    v = paras["freqs"] * math.sin(paras["dir"])
    # Computations.
    carrier = G2_carrier_real(x=x, y=y, freqx=u, freqy=v, phase=paras["phase"])
    envelope = G2_envelope(x=x, y=y, amp=paras["amp"], sdev=paras["sdev"])
    out = carrier * envelope
    return out

def G2_fil_imag(x, y, paras):
    '''
    Custom-built filter response (imaginary part).
    '''
    # Spatial frequency in polar coordinates.
    u = paras["freqs"] * math.cos(paras["dir"])
    v = paras["freqs"] * math.sin(paras["dir"])
    # Computations.
    carrier = G2_carrier_imag(x=x, y=y, freqx=u, freqy=v, phase=paras["phase"])
    envelope = G2_envelope(x=x, y=y, amp=paras["amp"], sdev=paras["sdev"])
    out = carrier * envelope
    return out

定義域が2次元なので、これらの関数のグラフも(色の濃淡などを使って)表すことができる。

In [11]:
myfig = plt.figure(figsize=(18,4))

y0 = math.floor(PIX_H/2)
x0 = math.floor(PIX_W/2)
y_inputs, x_inputs = np.mgrid[-y0:(y0+1), -x0:(x0+1)]

# Store pixel values (of envelope, carrier, filter) for plotting via imshow.
out_envelope = G2_envelope(x=x_inputs, y=y_inputs,
                             amp=myparas["amp"],
                             sdev=myparas["sdev"])
out_carrier = G2_carrier_real(x=x_inputs,
                              y=y_inputs,
                              freqx=myparas["freqs"]*math.cos(myparas["dir"]),
                              freqy=myparas["freqs"]*math.sin(myparas["dir"]),
                              phase=myparas["phase"])
out_filter = G2_fil_real(x=x_inputs, y=y_inputs, paras=myparas)


ax_carrier = myfig.add_subplot(1,3,1)
plt.title("Sinusoidal carrier")
topass = ax_carrier.imshow(out_carrier, cmap=plt.cm.BuPu_r)
plt.colorbar(topass)

ax_envelope = myfig.add_subplot(1,3,2)
plt.title("Gaussian envelope")
topass = ax_envelope.imshow(out_envelope, cmap=plt.cm.BuPu_r)
plt.colorbar(topass)

ax_filter = myfig.add_subplot(1,3,3)
plt.title("Gabor filter (2D)")
topass = ax_filter.imshow(out_filter, cmap=plt.cm.BuPu_r)
plt.colorbar(topass)

plt.show()

練習問題 (C):

  1. フィルタの「向き」を制御するパラメータdir($\theta$に対応)をゼロ、$\pi/4$、$\pi/2$など色々な値にしてみてください。数値を変えることで、フィルタの効能がどう変わると思われるか。

  2. 空間周波数パラメータfreqsを変えたときに、まず画像全体にわたって何周するか確認してみてください。踏まえて、「ヨコ」方向のみに効くようにdirを変え、ちょうど5周するようにfreqsも変えること。同様に、10周、25周するように。

  3. 直前のコードでは実部(carrier_realG2_fil_real)しか使っていない。これらを虚部に変えることで、フィルタの出力がどのように変わるか。

  4. 標準偏差のパラメータを調整し、フィルタの有効な半径が(視認できる範囲で)画像幅の半分になるようにすること。同様に、画像幅の1/4、1/8。


2Dフィルタの応用

基礎をおさえた上で、2次元版のガボールフィルタ(空間フィルタ)を応用していこう。画素値そのものが、2次元平面を定義域とする信号の値に相当する。

まずは簡単な応用例から入る。デジタル画像の読み書きを楽にしてくれるライブラリとして、imageio ( http://imageio.github.io/ )が著名である。画像をすでにローカルに保存してることを前提に、次の事例に使う画像を容易に読み込めることを確認しておく。

In [12]:
import imageio

# Read images from file.
im_cat = imageio.imread("img/chelsea.png")
im_boy = imageio.imread("img/bishop.png")
im_parrots = imageio.imread("img/parrots.png")

# Inspect (and plot) the images to ensure they have been read as we expect.
print("Shape:", im_cat.shape, "Type:", type(im_cat), "First value (r,g,b):", im_cat[0,0,:])
print("Shape:", im_boy.shape, "Type:", type(im_boy), "First value (r,g,b):", im_boy[0,0,:])
print("Shape:", im_parrots.shape, "Type:", type(im_parrots), "First value (r,g,b):", im_parrots[0,0,:])

myfig = plt.figure(figsize=(18,4))
ax_cat = myfig.add_subplot(1,3,1)
ax_boy = myfig.add_subplot(1,3,2)
ax_parrots = myfig.add_subplot(1,3,3)
ax_cat.imshow(im_cat)
ax_boy.imshow(im_boy)
ax_parrots.imshow(im_parrots)

plt.show()
Shape: (300, 451, 3) Type: <class 'imageio.core.util.Image'> First value (r,g,b): [143 120 104]
Shape: (128, 128, 3) Type: <class 'imageio.core.util.Image'> First value (r,g,b): [180 116  33]
Shape: (633, 801, 3) Type: <class 'imageio.core.util.Image'> First value (r,g,b): [168 165 160]

注目すべきは、imageioが自前の画像オブジェクトを使っていること(ndarray型ではない)。その利点は、RGBチャネルがきちんと認識され、カラー画像として表示されることである。

「色」を計算機で表現することは決して自明ではないが、出発点として普通のPNG形式画像のRGBチャネルを取り出していく。

In [13]:
im = imageio.imread("img/chelsea.png")

myfig = plt.figure(figsize=(18,4))
ax_R = myfig.add_subplot(1,3,1)
plt.title("Red channel")
ax_G = myfig.add_subplot(1,3,2)
plt.title("Green channel")
ax_B = myfig.add_subplot(1,3,3)
plt.title("Blue channel")

ax_R.imshow(im[:,:,0], plt.get_cmap('gray'))
ax_G.imshow(im[:,:,1], cmap=plt.get_cmap('gray'))
ax_B.imshow(im[:,:,2], cmap=plt.get_cmap('gray'))
plt.show()

各チャネルの度合い(白いほど高い)はこれで十分わかるが、やや強引にそれらが相当する色に結びつけようと思うと、カラーマップは自由に帰ることができる。

In [14]:
myfig = plt.figure(figsize=(18,4))
ax_R = myfig.add_subplot(1,3,1)
plt.title("Red channel")
ax_G = myfig.add_subplot(1,3,2)
plt.title("Green channel")
ax_B = myfig.add_subplot(1,3,3)
plt.title("Blue channel")

ax_R.imshow(im[:,:,0], cmap=plt.cm.Reds)
ax_G.imshow(im[:,:,1], cmap=plt.cm.Greens)
ax_B.imshow(im[:,:,2], cmap=plt.cm.Blues)
plt.show()