- 講師：山崎俊彦
- 参考書：CG-ARTS協会発行「ディジタル画像処理」
- 参考書：R. Szeliski, Computer Vision Algorithms and Applications, Springer (PDF版はインターネット上で無料公開)
領域の特徴量
In [1]:
import cv2
import numpy as np
import matplotlib.pyplot as plt
（２次元）フーリエ変換
- 画像は、フーリエ変換によって周波数領域における別の形で表現される。
- フーリエ変換の結果は、画像に含まれるそれぞれの周波数成分を表している。
- フーリエ変換後の各周波数成分の大きさを成分ごとに変えることにより、元の画像の性質を変えることができる。（周波数フィルタリング(frequency filtering)）
- ローパスフィルタ(lowpass fileter)：画像に含まれる空間周波数成分のうち、低周波数成分は残し、高周波数成分は除去するようなフィルタ。
- ハイパスフィルタ(highpass filter)：上の逆。高周波数成分は残し、低周波数成分を除去するようなフィルタ。
- バンドパスフィルタ(bandpass filter)：ある中間的な周波数の範囲を残すようなフィルタ。
- 高域強調フィルタ(high-emphasis filter)：画像の低周波数成分はそのまま保ちつつ、高周波数成分を強調するフィルタ。
$1$ 次元フーリエ変換¶
簡単のため、まずは一次元フーリエ変換をみる。
| 名前 | 数式 | 説明 |
|---|---|---|
| フーリエ変換 | $$F(\omega) = \int_{-\infty}^{\infty}f(t)e^{-i\omega t}dt$$ | 数直線上の表現（空間領域）から周波数による表現（周波数領域）への変換 |
| 逆フーリエ変換 | $$f(x) = \int_{-\infty}^{\infty}F(\omega)e^{i\omega x}d\omega$$ | 周波数領域にある関数を元の空間領域に戻す変換 |
実装¶
二つの正弦波を重ねた波から、それぞれの成分を取り出すことができるか確かめる。なお、ここでは numpy.fft.fft を利用しています。
これは、高速フーリエ変換(Fast Fourier Transform)と呼ばれる、高速に「離散フーリエ変換」を行うアルゴリズムを実装したものです。
In [2]:
N = 10000
dt = 0.01
time = np.arange(0, N*dt, dt)
freq = np.linspace(0, 1.0/dt, N)
In [3]:
def sin(A,f,p):
"""
@param A: Amplitude
@param f: Frequency[Hz]
@param p: phase
"""
func = lambda t:A*np.sin(2*np.pi*f*t + p1)
return func
In [4]:
# parameters
f1, f2 = 5, 8 # frequency[Hz]
A1, A2 = 5, 3 # Amplitude
p1, p2 = 0, 1/2*np.pi # phase
In [5]:
y1 = sin(A1,f1,p1)(time)
y2 = sin(A2,f2,p2)(time)
y = y1 + y2 # superposition
In [6]:
# Computes the discrete Fourier transform (DFT) of a sequence
yf = np.fft.fft(y)/(N/2) # Fast Fourier transform.
In [7]:
fig = plt.figure(figsize=(8,4))
ax1=fig.add_subplot(2,1,1)
ax1.set_ylabel("amplitude")
ax1.set_xlabel("time")
ax1.set_xlim(0,1)
ax1.plot(time,y, label="y(y1+y2)", color="blue")
ax1.plot(time,y1, label=f"y1(A={A1},f={f1})", alpha=0.3, color="green")
ax1.plot(time,y2, label=f"y2(A={A2},f={f2})", alpha=0.3, color="red")
ax1.axhline(0, alpha=0.3, color="black")
ax1.legend()
ax2=fig.add_subplot(2,1,2)
ax2.set_ylabel("amplitude")
ax2.set_xlabel("frequency")
ax2.set_xlim(0,10)
ax2.plot(freq, abs(yf), color="blue")
plt.tight_layout()
plt.show()
見事に、二つの周波数成分を取り出すことができました。
$2$ 次元フーリエ変換¶
※ 一般に画像は2次元空間で定義されたある関数 $f(x,y)$ として表現することができる！
| 名前 | 数式 |
|---|---|
| フーリエ変換 | $$F(u,v) = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x,y)\exp\left\{-j2\pi\left(ux+vy\right)\right\}dxdy$$ |
| 逆フーリエ変換 | $$f(x,y) = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}F(u,v)\exp\left\{j2\pi\left(ux+vy\right)\right\}dudv$$ |
ここで、以下の値に注目することがよくある。
| 名前 | 数式 | 説明 |
|---|---|---|
| 振幅スペクトル(amplitude spectrum) | $$\mid F(u,v)\mid = \sqrt{\mathrm{Re}\left\{F(u,v)\right\}^2 + \mathrm{Im}\left\{F(u,v)\right\}^2}$$ | 絶対値 |
| 位相スペクトル(phase spectrum) | $$\arg\left\{F(u,v)\right\} = \tan^{-1}\frac{\mathrm{Im}\left\{F(u,v)\right\}}{\mathrm{Re}\left\{F(u,v)\right\}}$$ | 偏角 |
| パワースペクトル(power spectrum) | $$\mid F(u,v)\mid^2$$ |
In [8]:
img = cv2.imread('lenna.png',0)
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
magnitude_spectrum = 20*np.log(np.abs(fshift))
plt.subplot(121), plt.imshow(img, cmap='gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.xlabel('x'), plt.ylabel('y')
plt.subplot(122), plt.imshow(magnitude_spectrum, cmap='gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.xlabel('u'), plt.ylabel('v')
plt.show()
ガボールフィルタ(Gabor Filter)
- 正弦波とガウス関数の積で表される。（例：2次元ガボールフィルタ）
$$g(x,y,\lambda,\phi) = \exp\left(-\frac{x^2+y^2}{2\sigma^2}\right)\exp\left(2\pi\lambda i\left(x\cos\phi + y\sin\phi\right)\right)$$
- 実数部： $$\mathrm{Re}\left\{g(x,y,\lambda,\phi)\right\} = \exp\left(-\frac{x^2+y^2}{2\sigma^2}\right)\cos\left(2\pi\lambda i\left(x\cos\phi + y\sin\phi\right)\right)$$
- 虚数部： $$\mathrm{Im}\left\{g(x,y,\lambda,\phi)\right\} = \exp\left(-\frac{x^2+y^2}{2\sigma^2}\right)\sin\left(2\pi\lambda i\left(x\cos\phi + y\sin\phi\right)\right)$$
- フーリエ変換では、所望の周波数成分が領域のどこにあるのかわからないのに対し、ガボールフィルタは、「領域のどの位置にどの方向の周波数が、どの強度であるか」を知ることができる。
- パラメータを操作して様々な方向と強度を持つフィルタ（フィルタバンク(filter bank)）を多数準備し、特徴抽出を行う。また、ここで抽出された特徴量群はtextonと呼ばれる。
- 入力された画像に対してガボールフィルタを畳み込んで（たたみ込み＝二次元の内積＝類似度を反映）、特徴量 $G(x,y,\lambda,\phi)$ を求める式は以下となる。 $$G(x,y,\lambda,\phi) = \sum_u\sum_vI(x,y)\ast g(x-u, y-v, \lambda, \phi)$$
OpenCVの cv2.getGaborKerne を用いることで簡単にガボールフィルタを生成することができる、
| 引数 | 説明 |
|---|---|
ksize |
Size of the filter returned. |
sigma |
Standard deviation of the gaussian envelope. |
theta |
Orientation of the normal to the parallel stripes of a Gabor function. |
lambd |
Wavelength of the sinusoidal factor. |
gamma |
Spatial aspect ratio. |
psi |
Phase offset. |
In [9]:
default_gabor = {
"ksize": (50,50),
"sigma": 4,
"theta": 0,
"lambd": 10,
"gamma": 0.5,
"psi": 0,
}
In [10]:
gabor = cv2.getGaborKernel(**default_gabor)
plt.imshow(gabor, cmap="gray")
plt.title("default")
plt.show()
In [11]:
def argsGaborRelation(key, vals):
global default_gabor
defaul_val = default_gabor[key]
N = len(vals)
n_col = 4
n_row = N//n_col if N%n_col==0 else N//n_col+1
fig = plt.figure(figsize=(4*n_col,4*n_row))
for i,val in enumerate(vals):
ax=fig.add_subplot(n_row,n_col,i+1)
default_gabor[key] = val
gabor = cv2.getGaborKernel(**default_gabor)
ax.imshow(gabor, cmap="gray")
ax.set_title(f"{key}={val}")
ax.set_xticks([]), ax.set_yticks([])
plt.show()
default_gabor[key] = defaul_val
In [12]:
argsGaborRelation(key="ksize", vals=[(30,30),(20,40),(40,20),(50,50)])
In [13]:
argsGaborRelation(key="sigma", vals=[x for x in range(8)])
In [14]:
argsGaborRelation(key="theta", vals=[np.radians(x) for x in range(0,180,45)])
In [15]:
argsGaborRelation(key="gamma", vals=[x for x in range(4)])
In [16]:
argsGaborRelation(key="psi", vals=[np.radians(x) for x in range(0,360,45)])
In [17]:
img = cv2.imread('lenna.png',0)
color_img = cv2.imread('lenna.png')
In [18]:
def argsGaborSuperposition(key, vals):
global default_gabor, img
defaul_val = default_gabor[key]
N = len(vals)*2
n_col = 8
n_row = N//n_col if N%n_col==0 else N//n_col+1
fig = plt.figure(figsize=(4*n_col,4*n_row))
dests=0
for i,val in enumerate(vals):
default_gabor[key] = val
gabor = cv2.getGaborKernel(**default_gabor)
dest = cv2.filter2D(img, cv2.CV_32F, gabor)
dests+=dest
ax_dst=fig.add_subplot(n_row,n_col,2*i+1)
ax_dst.imshow(np.power(dest,2), cmap="gray"),ax_dst.set_xticks([]), ax_dst.set_yticks([])
ax_gab=fig.add_subplot(n_row,n_col,2*i+2)
ax_gab.imshow(gabor, cmap="gray"),ax_gab.set_xticks([]), ax_gab.set_yticks([])
# Superposition
fig = plt.figure(figsize=(8,16))
ax_org=fig.add_subplot(1,2,1)
ax_org.imshow(cv2.cvtColor(color_img, cv2.COLOR_BGR2RGB))
ax_org.set_title("original image"),ax_org.set_xticks([]), ax_org.set_yticks([])
ax_dst=fig.add_subplot(1,2,2)
ax_dst.imshow(dests, cmap="gray")
ax_dst.set_title(f"Superposition of dst (key={key})"),ax_dst.set_xticks([]), ax_dst.set_yticks([])
plt.show()
default_gabor[key] = defaul_val
In [19]:
argsGaborSuperposition(key="theta", vals=[np.radians(x) for x in np.linspace(0,180,8)])
※ 入力画像の８分割した全ての方向のエッジを重ね合わせると、元の画像をかなり高い精度で復元していることがわかる。
同時生起確率
- テクスチャの統計的特徴量を求める。
- 離れた２つの場所にある画素対の値から、画素値の一様性・方向性・コントラストなどの性質を表す特徴量を求める。
- ある画素 $i$ と、$i$ から離れた位置にある画素 $j$ の画素対を考え、画素 $i$ と画素 $j$ の相対的な位置を $\delta=\langle d,\theta\rangle$ とする。 それぞれの画素値を $L_i,L_j$ とし、画素値の対 $(L_i,L_j)$ が生じる出現頻度で、ある同時生起行列 $\mathbf{H}_{\delta}(L_i,L_j)$ を考える。ここで、出現頻度の総数で $\mathbf{H}_{\delta}(L_1,L_2)$ を正規化し、確率に変換した同時生起行列を $\mathbf{P}_{\delta}(L_i,L_j)$ とする。
- 確率に変換した同時生起行列を用いて計算する主な特徴量には次の $7$ 種類がある。（ここで、画素値のレベル数は $L$ である。）他にも様々な統計量を計算できるが、テクスチャの物理的な意味と直接には結びつかない数値もある。
$$
\begin{aligned}
\mathrm{ENR} &= \sum_{L_i=0}^{L-1}\sum_{L_j=0}^{L-1}\left\{P_{\delta}(L_i,L_j)\right\}^2 & (\text{energy})\\
\mathrm{CNT} &= \sum_{k=0}^{L-1}k^2P_{x-y}(k) & (\text{contrast})\\
\mathrm{CRR} &= \frac{\sum_{L_i=0}^{L-1}\sum_{L_j=0}^{L-1}L_iL_jP_{\delta}(L_i,L_j) - \mu_x\mu_y}{\delta_x\delta_y} & (\text{correlation})\\
\mathrm{VAR} &= \sum_{L_i=0}^{L-1}\sum_{L_j=0}^{L-1}\left(L_i-\mu_x\right)^2P_{\delta}(L_i,L_j) & (\text{variance})\\
\mathrm{EPY} &= -\sum_{L_i=0}^{L-1}\sum_{L_j=0}^{L-1}P_{\delta}(L_i,L_j)\log\left\{P_{\delta}(L_i,L_j)\right\} & (\text{entropy})\\
\mathrm{SEP} &= -\sum_{k=0}^{2L-2}P_{x+y}(k)\log\left\{P_{x+y}(k)\right\} & (\text{sum entropy})\\
\mathrm{IDM} &= \sum_{L_i=0}^{L-1}\sum_{L_j=0}^{L-1}\frac{1}{1+(L_i-L_j)^2}P_{\delta}(L_i,L_j) & (\text{inverse difference moment})\\
\end{aligned}
$$
ただし、
$$ \mu_x = \sum_{L_i=0}^{L-1}L_iP_x(L_i), \quad \mu_y = \sum_{L_j=0}^{L-1}L_jP_y(L_j)\\ \delta_x^2 = \sum_{L_i=0}^{L-1}(L_i-\mu_x)^2P_x(L_i), \quad \delta_y^2 = \sum_{L_j=0}^{L-1}(L_j-\mu_y)^2P_y(L_j)\\ P_x(L_i) = \sum_{L_j=0}^{L-1}P_{\delta}(L_i,L_j), \quad P_y(L_j) = \sum_{L_i=0}^{L-1}P_{\delta}(L_i,L_j)\\ P_{x-y}(k) = \sum_{L_i=0}^{L-1}\sum_{L_j=0}^{L-1}P_{\delta}(L_i,L_j)\quad |L_i-L_j|= k\\ P_{x+y}(k) = \sum_{L_i=0}^{L-1}\sum_{L_j=0}^{L-1}P_{\delta}(L_i,L_j)\quad |L_i+L_j|= k\\ $$