本記事のゴール
スケールの変化や回転に強い画像の特徴を記述する SIFT (Scale Invariant Feature Transform) を支える Scale-space を理解し実装できるようになる。
関連して得られる知識
- Gaussian convolution による画像の平滑化
- Fourier 変換とサンプリング定理
使用する環境
| ソフトウェア名 | バージョン |
|---|---|
| Python | 3.4 or 3.5 |
| Pillow | 3.1.0 |
| Numpy | 1.10 |
| Scipy | 0.16.1 |
| Matplotlib | 1.5.0 |
注釈
本記事は Morning Project Samurai 第41回ミーティングの資料として書かれたものである。
お願い
これはドラフトバージョンです。ミスなどを発見された場合は、ご連絡いただけると幸いです。
Scale-space
定義
- 与えられた画像からスケールが異なる複数の画像を生成する。
- それら画像で構成される空間を Scale-space と呼ぶ。
モチベーション
任意の与えられる画像からある特定の物体の検出を考える。例えば市街地道路の動画の一場面 （画像）から人を検出することを考える。この時、人の着ている服の模様の詳細がわかるようなスケールの画像は必要ない。人の形はわかるが着ている服はどんなものかわからない程度のスケールの画像の方が良い。詳細すぎる情報はむしろノイズとなる。
上述のように、目的の物体を検出するのに適した画像のスケールが存在する。しかし、現実では与えられる画像がいつも適切なスケールとは限らない。そこで、与えられる画像からスケールが異なる複数の画像を生成して集めた Scale-space を構成し、その中から適切なスケールの画像を見つけ物体検出を行う。
スケールとスケールアップ
Scale-space においてスケールアップとは、Gaussian convolution と呼ばれる手法を用いて元画像から詳細な情報を取り除くことである。スケールアップの度合いは Gaussian convolution のパラメタ $\sigma$ によって決定される。スケール $\sigma$ が大きくなるほどより多くの情報が元画像から取り除かれる。
例えば下図であれば、一番左のオリジナル画像は顔の細部や帽子のファーの毛まではっきりと見て取ることができる。一番右のスケール $\sigma = 3.2$ の画像では、ぼんやりと女性がおり目鼻立がはっきりしているなということはわかるが、それ以上のことはわからない。
構成方法
- 元画像（グレースケール）を準備。
- スケール $\sigma = k^{n}\sigma_{0}$ とし、$n$ の初期値を $0$ とする。
- 元画像を Gaussian convolution を用いてスケールアップ。
- $k^{n}=2.0$ となるまで $n$ を $1, 2, 3,...$ と増加させながら手順3を繰り返す。$k$ は $k^{n}=2.0$ となるまでに何枚の画像を作りたいかによって決定される。
- スケール $\sigma = 2.0\sigma_{0}$ の画像を縦と横が半分のサイズになるようダウンサンプリング。
- ダウンサンプリングした画像を元画像とし、$n=1$ として手順 3 ヘ戻る。
詳しい構成方法とそのプログラムは、 Gaussian convolution や Fourier 変換、サンプリング定理などを一通り学習したのち述べる。
Gaussian convolution
- Scale-space では元画像からスケール $\sigma$ の画像を生成する際に用いられる。
- 一般的な画像処理の文脈では画像のノイズ除去やぼかしに用いられる。
- 数学的には Gaussian kernel を用いて任意の関数の重みつき平均の系列を求める。
1次元 Gaussian convolution
1次元 Gaussian convolution は次式で表される。
F(x, \sigma)=\int_{-\infty}^{\infty} f(u) \frac{1}{\sqrt{2 \pi \sigma^{2}}}e^{-\frac{(u-x)^2}{2\sigma^{2}}} du
$f$ は元信号である。$\frac{1}{\sqrt{2 \pi \sigma^{2}}}e^{-\frac{(u-x)^2}{2\sigma^{2}}}$ は Gaussian kernel と呼ばれる。$\sigma$ は Gaussian kernel の形状を決定するパラメタである。 Gaussian kernel は中心 $x$ で最大値を取る釣鐘のような形をしており、その高さと広がり具合が $\sigma$ で決定される。
Gaussian kernel
試しに numpy と scipy 、matplotlib を用いて Gaussian kernel を描画してみよう。
import numpy as np
from scipy.signal import gaussian
from matplotlib import pyplot as plt
if __name__ == '__main__':
npoints = 101
sigmas = np.array([6.0, 12.0])
for i, sigma in enumerate(sigmas):
y = gaussian(npoints, sigma) / (np.sqrt(2.0 * np.pi) * sigma)
plt.subplot(len(sigmas), 1, i + 1)
plt.title('sigma = %s' % sigma)
plt.ylim(ymax=0.08)
plt.plot(np.arange(npoints/2 - npoints, npoints/2, dtype=np.int), y)
plt.tight_layout()
plt.show()
上記のプログラムを実行すれば、下図が結果として得られる。
[Gaussian kernel の性質を後日追記]
解釈と適用例
Gaussian convolution の式は「元信号 $f$ に中心 $x$ 広がり $\sigma$ の重み (Gaussian kernel) を掛けてその総和をとったもの（加重平均）を $F(x, \sigma)$ で表す」と解釈できる。これより、$F$ は「元信号 $f(u)$ の $u=x$ 近傍の情報を $\sigma$ で決定づけられる強さで色濃く反映しつつ $f$ 全体の情報を持った点 $F(x, \sigma)$ で構成される新たな信号」であると解釈できる。これは「$F$ は元信号 $f$ を強度 $\sigma$ で平滑化した（ぼかした、ノイズ除去した）信号」だということである。Scale-space の文脈では「 $F$ は $f$ をスケール $\sigma$ にスケールアップした（詳細な情報を取り除いた）信号」だと言える。
Gaussian kernel は前節で見たように $\sigma$ が大きくなるほど平坦になる。これは、「$\sigma$ が大きくなるほど $F(x, \sigma)$ に含まれる元信号 $f$ の $u=x$ の近傍の情報の色が薄れる」ことを意味し、$\sigma$ をどんどん大きくしていけば最終的に $F$ は全ての $x$ において同じような値をとるようになる。Scale-space の文脈では「スケール $\sigma$ を大きくするほど巨視的な視点で信号を眺め（詳細な情報は失われ）、最終的に信号に含まれる情報をすべて同一視するようになる」と解釈できる。
ここまでの話をプログラムで可視化してみよう。
import numpy as np
from scipy.ndimage.filters import gaussian_filter1d
from matplotlib import pyplot as plt
if __name__ == '__main__':
x = np.arange(0, 100)
f = np.sin(0.5 * x) + np.random.normal(0, 6.0, x.shape)
sigmas = [1.6, 3.2, 6.4]
plt.subplot(len(sigmas) + 1, 1, 1)
plt.ylim(ymin=np.min(f), ymax=np.max(f))
plt.title('Original')
plt.plot(f)
for i, sigma in enumerate(sigmas):
plt.subplot(len(sigmas) + 1, 1, 2 + i)
plt.title('Sigma = %s' % sigma)
plt.plot(gaussian_filter1d(f, sigma))
plt.tight_layout()
plt.show()
上記プログラムを実行すると下図が得られる。
$\sigma$ が大きくなる（スケールアップする）につれ、$F$ から信号 $f$ の含む細かな情報がどんどん失われていく様がわかる。
2次元 Gaussian convolution
画像は2次元信号である。そこで、Gaussian convolution を用いて画像をスケールアップするためには2次元信号に対する Gaussian convolution が必要となる。以降、単に Gaussian convolution と書く時は、2次元 Gaussian convolution を指すこととする。
2次元信号に対する Gaussian convolution は次式で表される。
F(x, y, \sigma) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(u, v) \frac{1}{2\sigma^{2}} e^{-\frac{(u - x)^2 + (v-y)^{2}}{\sigma^{2}}} dudv
$\frac{1}{2\sigma^{2}} e^{-\frac{(u - x)^2 + (v-y)^{2}}{\sigma^{2}}}$ が Gaussian kernel である。この Gaussian kernel は立体の釣鐘型であり、$(x, y)$ で最大値を取り $\sigma$ はその広がりを表す。
Gaussian kernel
プログラムで Gaussian kernel を描画してみよう。
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def gaussian_kernel_2d(x, y, sigma):
return np.exp(-(np.power(x/sigma, 2) + np.power(y/sigma, 2)))/(2 * np.power(sigma, 2))
if __name__ == '__main__':
xrange = np.arange(-10.0, 10.5, 0.5)
yrange = np.arange(-10.0, 10.5, 0.5)
kernel_values = np.zeros(shape=(len(yrange), len(xrange)))
sigmas = [1.6, 2.4]
fig = plt.figure()
X, Y = np.meshgrid(xrange, yrange)
for i, sigma in enumerate(sigmas):
ax = fig.add_subplot(len(sigmas), 1, 1 + i, projection='3d')
plt.title('Sigma = %s' % sigma)
ax.set_zlim(zmax=0.2)
ax.plot_wireframe(X, Y, gaussian_kernel_2d(X, Y, sigma))
plt.show()
上記プログラムを実行すると下図を得られる。マウスで視点を変えることが可能である。
離散化 Gaussian convolution
Gaussian filter
コンピュータは離散的な値しか扱うことができず、またその範囲は有限である。Gaussian kernel は連続な関数で定義域は $\pm \infty$ である。そこで、コンピュータで扱える形にするため、Gaussian kernel を有限個の数値列で近似して離散化する必要がある。
$s$ 行 $t$ 列の格子を考える。以降この格子を $s \times t$ の filter という。これはプログラム上で下記のような2次元配列で表現される。
filter = np.zeros(shape=(s, t))
変数 filter のすべての要素に $u=0$ の時の Gaussian kernel をよく表現するよう値を設定できれば、$u=0$ の Gaussian kernel を離散化した Gaussian filter $g$ が作成できる。ここでは、$g$ の $k$ 行 $l$ 列の要素 $g(k,l)$ = $\frac{1}{\alpha} e^{-\frac{l^{2} + k^{2}}{\sigma^{2}}}$ と設定する。ここで $\alpha$ は $g$ の全ての要素を足した時にその合計が $1$ となるように正規化する定数である。
Gaussian filter を用いた Gaussian convolution
次式は Gaussian filter を用いた Gaussian convolution である。これを離散化 Gaussian convolution と呼ぶ。以降、単に Gaussian convolution というときは、この離散化 Gaussian convolution を指すものとする。
F(i, j) = \sum_{k=0}^{s}\sum_{l=0}^{t} f(i+k-[\frac{s}{2}], j+l-[\frac{t}{2}]) g(k, l)
- $f$: $m \times n$ の元画像
- $g$: $s \times t$ の Gaussian filter
- $F$: $m \times n$ の出力画像
上式の通りに素直にプログラムを作成すると下記のようになる。
for k in range(s):
for l in range(t):
F[i, j] = f[i + k - int(s/2), j + l - int(t/2)] * g(k, l)
この場合出力画像 $F$ の1画素を求めるために $s \times t$ 回の演算が必要であり、出力画像全ての画素を求めるためには $m \times n \times s \times t$ 回の演算が必要となる。$320 \times 240$ の画像と $5 \times 5$ の Gaussian filter で convolution した場合、合計で $1,920,000$ 回の演算が必要となる。この演算回数がクリティカルであるか否かはマシンスペックやアプリケーションによるが、私の MacBook Air (13-inch, Early 2015) 上ではクリティカルである。
高速化
Gaussian convolution は、$s \times t$ の Gaussian filter を用いた Gaussian convolution の出力画像 $F$ が、次の 2 ステップで出力される画像 $F_{1}$ で表されるという性質を持っている。
- $1 \times t$ の 1 次元 Gaussian filter を元画像 $f$ の各行と convolution して画像 $F_{0}$ を出力。
- $s \times 1$ の 1 次元 Gaussian filter を $F_{0}$ の各列と convolution して画像 $F_{1}$ を出力。
これは次のように証明される。
\begin{align}
F(i, j) &=\sum_{k=0}^{s}\sum_{l=0}^{t} f(i+k-[\frac{s}{2}], j+l-[\frac{t}{2}]) g(k, l)\\
&= \sum_{k=0}^{s}\sum_{l=0}^{t} f(i+k-[\frac{s}{2}], j+l-[\frac{t}{2}]) \frac{1}{\alpha}e^{-\frac{l^{2} + k^{2}}{\sigma^{2}}}\\
&= \sum_{k=0}^{s}\sum_{l=0}^{t} f(i+k-[\frac{s}{2}], j+l-[\frac{t}{2}]) \frac{1}{\alpha}e^{-\frac{l^{2}}{\sigma^{2}}}e^{-\frac{k^{2}}{\sigma^{2}}}\\
&=\frac{\alpha_{0}\alpha_{1}}{\alpha} \sum_{k=0}^{s}\frac{1}{\alpha_{1}}e^{\frac{-k^{2}}{\sigma^{2}}}\sum_{l=0}^{t} f(i+k-[\frac{s}{2}], j+l-[\frac{t}{2}]) \frac{1}{\alpha_{0}}e^{-\frac{l^{2}}{\sigma^{2}}}\\
&=\frac{\alpha_{0}\alpha_{1}}{\alpha} \sum_{k=0}^{s}\frac{1}{\alpha_{1}}e^{\frac{-k^{2}}{\sigma^{2}}}\sum_{l=0}^{t} f(i+k-[\frac{s}{2}], j+l-[\frac{t}{2}]) g_{0}(l)\\
&=\frac{\alpha_{0}\alpha_{1}}{\alpha} \sum_{k=0}^{s}\frac{1}{\alpha_{1}}e^{-\frac{k^{2}}{\sigma^{2}}} F_{0}(i, j)\\
&=\frac{\alpha_{0}\alpha_{1}}{\alpha} \sum_{k=0}^{s} F_{0}(i, j) g_{1}(k)\\
&=\frac{\alpha_{0}\alpha_{1}}{\alpha} F_{1}(i, j)
\end{align}
- $g_{0}$: $1 \times t$ の Gaussian filter
- $g_{1}$: $s \times 1$ の Gaussian filter
- $\alpha$: $g$ の正規化定数
- $\alpha_{0}$: $g_{0}$ の正規化定数
- $\alpha$: $g_{1}$ の正規化定数
このようにすることで演算回数は $s \times m + t \times n$ となる。$320 \times 240$ の画像と $5 \times 5$ の Gaussian filter の convolution の場合、演算回数は $768,000$ 回となる。高速化により演算回数を $60%$ 減らすことができる。
画像への適用
実際に画像に Gaussian convolution を適用してみよう。今回は Pillow, numpy, scipy を用いる。下記プログラムではスケール $\sigma = 1.6$ の画像とスケール $\sigma = 3.2$ の画像を生成している。
from PIL import Image
import numpy as np
from scipy.ndimage.filters import gaussian_filter
from matplotlib import pyplot as plt
if __name__ == '__main__':
orig_image = Image.open('lena.jpg').convert('L')
orig_image = np.array(orig_image, dtype=np.uint8)
sigmas = [1.6, 3.2]
plt.subplot(1, len(sigmas) + 1, 1)
plt.title('Orig')
plt.imshow(orig_image, cmap='Greys_r')
for i, sigma in enumerate(sigmas):
plt.subplot(1, len(sigmas) + 1, 2 + i)
plt.title('Sigma=%s' % sigma)
plt.imshow(gaussian_filter(orig_image, sigma), cmap='Greys_r')
plt.tight_layout()
plt.savefig('gausian_convolution_2d_examples.png')
plt.show()
上記プログラムを実行すると下図が得られる。
画像を扱う有名なライブラリに OpenCV がある。それがインストールされている人は、それを使っても良い。また、Gaussian convolution を自力でインプリメントしてみるのも良い。
ここまでで元画像に対して任意のスケール $\sigma$ の画像を生成できるようになった。
Scale-space 1st octave の構築
いよいよ Scale-space の構築を行う。章題中の 1st octave ついては後述する。
構築の流れ
- ベースとなるスケール $\sigma_{0}$ を決定。
- スケールが $k^{n}\sigma$ の画像を $n=0,...,s$ について生成する。$s$ は scale-space 1st octave の分割数。 $k=2^{1/s}$ で表される。
Scale-space 1st octave を構築するプログラム
from PIL import Image
import numpy as np
from scipy.ndimage.filters import gaussian_filter
from matplotlib import pyplot as plt
if __name__ == '__main__':
orig_image = Image.open('lena.jpg').convert('L')
orig_image = np.array(orig_image, dtype=np.uint8)
sigma = 1.6
s = 3
k = np.power(2, 1/s)
scale_space = []
for n in range(s + 1):
scale_space.append(gaussian_filter(orig_image, np.power(k, n) * sigma))
for n, img in enumerate(scale_space):
plt.subplot(1, len(scale_space), 1 + n)
plt.title('Sigma=%s' % np.round(np.power(k, n) * sigma, 2))
plt.imshow(img, cmap='Greys_r')
plt.tight_layout()
plt.show()
上記プログラムを実行すると変数 scale_space にスケールが $1.6$ から $3.2$ までの4枚の画像が格納され、下図が表示される。
この $k\sigma$ から $2k\sigma$ までの画像で構成されるスペースを 1st octave と呼ぶ。
Fourier 変換とサンプリング定理
2nd octave 以降を作成するにあたって、Fourier 変換とサンプリング定理について触れておく必要がある。ただし、今回は深入りせずさらっと流すのみにとどめる。
Fourier 変換
Fourier 変換は画像をピクセル領域から周波数領域へ変換するものである。
1次元 Fourier 変換の式
F(j\omega) = \int_{-\infty}^{\infty} f(x) e^{-j\omega x} dx
- $f$: 場所領域の信号（例えば細長いスリットを通過する光の明度）
- $F$: $f$ の周波数領域での表現
- $\omega$: 角周波数 ($2\pi\omega_{0}$、ここで $\omega_{0}$ は周波数)
- $e^{-j\omega x} = \cos(\omega x) - j \sin(\omega x) $
- $j$: 虚数単位
1次元 Fourier 変換の適用例
numpy, scipy, matplotlib を用いて FFT (Fast Fourier Transform) による Fourier 変換を行ってみる。下記プログラムでは、サンプリングレートを20, 元信号の角周波数を $2\pi$ としている。
import numpy as np
from scipy.fftpack import fft
from matplotlib import pyplot as plt
if __name__ == '__main__':
sampling_rate = 20
sampling_interval = 1.0/sampling_rate
x = np.arange(0, 1, sampling_interval)
omega0 = 1.0
omega = 2.0 * np.pi * omega0
f = np.sin(omega * x)
F = fft(f)
plt.plot(np.arange(-len(f)/2, len(f)/2), np.abs(np.concatenate((F[len(f)/2:], F[:len(f)/2]))), 'bo-')
plt.show()
上記プログラムを実行すると下図を得る。
この図から、元信号 $f$ には周波数 $1$ の信号が含まれていることがわかる。幾つかの周期関数を重ね合わせたものを Fourier 変換してみると良い。
[後日もう少し詳しく書く]
1次元信号のサンプリング定理
サンプリング定理は、「元信号 $f$ に含まれる信号の最大周期が $W$ であるとき、元信号のサンプリング周期 $T$ が $T\leq \frac{1}{W}$ であれば、サンプリング後の信号 $f_{d}$ から元信号 $f$ が完全に復元可能である」ことを主張する。
[証明は後日追加]
1次元信号のサンプリング定理の実験
サンプリング周期を $1/20$ とする。そうするとサンプリング定理より、最大周波数が $10$ の信号までなら Fourier 変換でその周波数情報を全てうまく取り出せるはずである。
下記プログラムで確認してみよう。
import numpy as np
from scipy.fftpack import fft
from matplotlib import pyplot as plt
if __name__ == '__main__':
sampling_rate = 20
sampling_interval = 1.0/sampling_rate
x = np.arange(0, 1, sampling_interval)
omega0s = [9, 10, 11, 30]
for i, omega0 in enumerate(omega0s):
omega = 2.0 * np.pi * omega0
f = np.sin(omega * x)
F = fft(f)
plt.subplot(len(omega0s), 1, 1 + i)
plt.plot(np.arange(-len(f)/2, len(f)/2), np.abs(np.concatenate((F[len(f)/2:], F[:len(f)/2]))), 'bo-')
plt.title('Max frequency = %s' % omega0)
plt.xlabel('Frequency')
plt.ylabel('Amplitude spectrum')
plt.tight_layout()
plt.show()
実行すると下図が得られる。この結果だと最大周波数が $9$ の信号で周波数情報がうまく取り出せていることがわかる（後日詳細追記）。
2次元 Fourier 変換
[後日追記]
2nd octave 以降の構築
ダウンサンプリング
概要
[後日追記]
ダウンサンプリングの妥当性証明
[後日追記]
2nd octave 構築プログラム
[後日追記]