pythonとopencvを使って画像処理を勉強していきます。
前回
python+opencvで画像処理の勉強3 領域に基づく濃淡変換(空間フィルタリング)
空間フィルタリングについて学びました。
今回は画像のフーリエ変換について学んでいきます。
画像のフーリエ変換
2次元フーリエ変換
画像を$f(x,y)$とするとき、そのフーリエ変換$F(u,v)$は以下の式で表される。
$$F(u,v)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x,y)\exp{(-j2\pi(ux+vy))dxdy}$$
$u,v$はそれぞれx,y方向の空間周波数であり、フーリエ変換によりx,yで表される空間（空間領域）から、u,vで表される別の空間（周波数領域）の別の形で表現されます。
また、$F(u,v)$はフーリエ逆変換により、元の画像$f(x,y)$に戻すことができます。
$$f(u,v)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}F(u,v)\exp{(j2\pi(ux+vy))dudv}$$
$F(u,v)$は複素数で、絶対値と偏角はそれぞれ
$$
|F(u,v)|=\sqrt{Re{F(u,v)}^2+Im{F(u,v)}^2}
$$
$$
arg{F(u,v)}=\tan^{-1}\frac{Im{F(u,v)}}{Re{F(u,v)}}
$$
と求めることができます。
ここで、$Re{F(u,v)}$、$Im{F(u,v)}$はそれぞれ$F(u,v)$の実部と虚部を表しています。
$|F(u,v)|$は画像の振幅スペクトル、$arg{F(u,v)}$は位相スペクトルと呼ばれ、$|F(u,v)|^2$はパワースペクトルと呼ばれます。
画像のフーリエ変換
各画像のフーリエ変換した結果を確認してみます。
imgs = [img_gray, img1_gray, img2_gray]
fig, ax = plt.subplots(2, 3, figsize=(12, 8))
for i in range(3):
dft = cv2.dft(np.float32(imgs[i]), flags = cv2.DFT_COMPLEX_OUTPUT)
# ゼロ周波数の成分を中心に移動
dft_shift = np.fft.fftshift(dft)
# パワースペクトル
magnitude_spectrum = 20*np.log(cv2.magnitude(dft_shift[:,:,0], dft_shift[:,:,1]))
"""
numpyによるフーリエ変換
f = np.fft.fft2(imgs[i])
fshift = np.fft.fftshift(f)
magnitude_spectrum = 20*np.log(np.abs(fshift))
"""
ax[0][i].imshow(imgs[i], 'gray')
ax[0][i].set_xticks([])
ax[0][i].set_yticks([])
ax[1][i].imshow(magnitude_spectrum, 'gray')
ax[1][i].set_xticks([])
ax[1][i].set_yticks([])
原点は直流成分を表しています。また、元画像のエッジと直交する方向にスペクトルが延びているのがわかります。
周波数フィルタリング
周波数フィルタリング
フーリエ変換後の各周波数成分の大きさを各成分毎に変えることにより、元の画像の性質を変化させることができます。
このような処理を周波数フィルタリングと呼びます。
フーリエ変換後の画像を$F(u,v)$、フィルタリングの出力を$G(u,v)$とするとき、周波数フィルタリングは以下の式で表されます。
$$G(u,v)=F(u,v)H(u,v)$$
ここで、$H(u,v)$は周波数フィルタを表しています。
$G(u,v)$の逆フーリエ変換を行うことでフィルタリングした結果の画像を見ることができます。
空間フィルタリングと周波数フィルタリングの関係
2つの関数$f_{1}(x,y),f_{2}(x,y)$に対し、フーリエ変換は以下に示す性質を持ちます。
F\{f_{1}(x,y)*f_{2}(x,y)\}=F\{f_{1}(x,y)\}F\{f_{2}(x,y)\}
ここで、$F{}$はフーリエ変換を表し、*は以下で定義される畳み込み積分を表しています。
f_{1}(x,y)*f_{2}(x,y) = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f_{1}(\xi,\eta)f_{2}(x-\xi,y-\eta)d\xi d\eta
この式は、2つの関数の畳み込み積分の結果をフーリエ変換したものは、それぞれの関数のフーリエ変換の籍に等しくなるということを意味しています。
ここで、入力画像を$f(x,y)$、フィルタ関数を$h(x,y)$としたとき、$f(x,y)*h(x,y)$は画像の空間フィルタリングを表します。
周波数フィルタリングでは、フーリエ変換とフーリエ逆変換をするため、その分計算時間は余分にかかります。
一方で、空間フィルタリングでは畳み込み積分をするところを、周波数フィルタリングではフィルタ関数との積で済むため、この部分では空間フィルタリングの方が時間を必要とします。
一般に空間フィルタが小さいときは空間フィルタリングが有利であり、その大きさが大きくなるに従い周波数フィルタリングが有利になると言えます。
ローパスフィルタ、ハイパスフィルタ、バンドパスフィルタ
ここで、いくつかの関数を定義しておきます。
- 円形のマスク画像を作る関数
- ガウス分布型のマスク画像を作る画像
- フーリエ変換をした画像にマスク処理を施して逆フーリエ変換をする関数
def make_mask(img, R, inv=False):
"""円形のマスク画像を作ります"""
height = img.shape[0]
width = img.shape[1]
center_w = height//2
center_h = width//2
if inv:
n = 0
filter_matrix = np.ones([height, width])
else:
n = 1
filter_matrix = np.zeros([height, width])
for i in range(0, height):
for j in range(0, width):
if (i-center_w)*(i-center_w) + (j-center_h)*(j-center_h) < R*R:
filter_matrix[i][j] = n
return filter_matrix
def make_Gauss_mask(x):
"""ガウス分布形のマスク画像を作ります"""
X = np.matrix(x).T
x_mu = X - mu
w1 = np.matmul(inv_sig, x_mu)
w2 = -0.5 * np.matmul(x_mu.T, w1)
N = 1/(np.sqrt(((2*np.pi)**2) * np.abs(det))) * np.exp(w2)
return N[0,0]
def masked_fft(img, mask):
dft = cv2.dft(np.float32(img), flags = cv2.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)*mask
magnitude_spectrum = 20*np.log(cv2.magnitude(dft_shift[:,:,0],dft_shift[:,:,1]))
magnitude_spectrum[magnitude_spectrum==-np.inf]=0
f_ishift = np.fft.ifftshift(dft_shift)
img_back = cv2.idft(f_ishift)
img_back = cv2.magnitude(img_back[:,:,0],img_back[:,:,1])
dft2 = cv2.dft(np.float32(img_back),flags = cv2.DFT_COMPLEX_OUTPUT)
dft_shift2 = np.fft.fftshift(dft2)
magnitude_spectrum2 = 20*np.log(cv2.magnitude(dft_shift2[:,:,0],dft_shift2[:,:,1]))
return magnitude_spectrum, img_back, magnitude_spectrum2
後で使うガウス分布のマスク画像を作成しておきます。
x1_range = img_gray.shape[1] // 2
x2_range = img_gray.shape[0] // 2
x1 = np.arange(-x1_range, x1_range)
x2 = np.arange(-x2_range, x2_range)
# ガウス分布の設定
sigmas = [100, 900, 1600]
M=[0, 0]
mu = np.matrix(M).T
g_masks = []
# 3パターンのガウス分布のマスク画像を作成
for i in range(3):
S=[[sigmas[i],0],[0,sigmas[i]]]
sigma = np.matrix(S)
inv_sig = np.linalg.inv(sigma)
det = np.linalg.det(sigma)
g2 = np.array([[make_Gauss_mask([i, j]) for i in x1] for j in x2])
g2 *= 1/g2.max()
g2 = np.array([g2, g2]).transpose(1, 2, 0)
g_masks.append(g2)
ローパスフィルタ
低周波数成分は残し、高周波数成分は除去するようなフィルタをローパスフィルタと呼びます。
結果としては平滑化された画像が得られます。遮断周波数が小さくなるほどより強く平滑化されます。
結果を確認してみます。
masks = np.array([np.array([make_mask(img_gray, r),
make_mask(img_gray, r)]).transpose(1, 2, 0) for r in [30, 50, 150]])
fig, ax = plt.subplots(2, 3, figsize=(12, 8))
for i in range(3):
spectrum_img1, img_ifft, spectrum_img2 = masked_fft(img_gray, masks[i])
ax[0][i].imshow(spectrum_img1, 'gray')
ax[0][i].set_xticks([])
ax[0][i].set_yticks([])
ax[1][i].imshow(img_ifft, 'gray')
ax[1][i].set_xticks([])
ax[1][i].set_yticks([])
ガウス分布型のローパスフィルタ
ガウス分布に従うローパスフィルタを用いた場合の結果を確認します。
上の結果では不自然な縞模様が表れていましたが、ガウス分布型ではより滑らかで自然な平滑化となっているのがわかります。
fig, ax = plt.subplots(2, 3, figsize=(12, 8))
for i in range(3):
spectrum_img1, img_ifft, spectrum_img2 = masked_fft(img_gray, g_masks[i])
ax[0][i].imshow(spectrum_img1, 'gray', vmin=0, vmax=300)
ax[0][i].set_xticks([])
ax[0][i].set_yticks([])
ax[1][i].imshow(img_ifft, 'gray')
ax[1][i].set_xticks([])
ax[1][i].set_yticks([])
空間フィルタリングによる平滑化との関係
ガウシアンフィルタ
$$
h_{g}(x,y) = \frac{1}{2\pi\sigma^2}\exp{(-\frac{x^2+y^2}{2\sigma^2})}
$$
のフーリエ変換は
$$
H_{g}(u,v)=\exp{(-2\pi^2\sigma^2(u^2+v^2))}
$$
と計算され、ガウス分布のフーリエ変換はガウス分布となることがわかります。
空間領域におけるガウシアンフィルタとガウス分布型のローパスフィルタは本質的に同じものとなります。
さらに、$h_{g}(x,y)$と$H_{g}(u,v)$の幅は反比例の関係となります。
次に平均化フィルタについて考えます。
平均化フィルタは以下のように表されます。
$$
h_{ave}(x,y)=\frac{1}{w^2}rect(\frac{x}{w}, \frac{y}{w})
$$
ただし、
rect(x, y)= \begin{cases}
1 & (|x|\leq\frac{1}{2}かつ|y|\leq\frac{1}{2}) \\
0 & (それ以外)
\end{cases}
そして$h_{ave}(x,y)$のフーリエ変換は以下のように計算されます。
$$
H_{ave}(u,v)=\frac{\sin{\pi wu}}{\pi wu}\frac{\sin{\pi wv}}{\pi wv}
$$
振幅スペクトルは周波数が高くなるにつれて小さくなるため、ローパスフィルタの性質を持つと言える。
ハイパスフィルタ
高周波数成分は残し、低周波数成分は除去するようなフィルタをハイパスフィルタと呼びます。
ハイパスフィルタは以下の式より、ローパスフィルタから作ることもできます。
$$H_{high}(u,v)=1-H_{low}(u,v)$$
masks_i = np.array([np.array([make_mask(img_gray, r, True), make_mask(img_gray, r, True)]).transpose(1, 2, 0) for r in [30, 50, 160]])
fig, ax = plt.subplots(2, 3, figsize=(12, 8))
for i in range(3):
spectrum_img1, img_ifft, spectrum_img2 = masked_fft(img_gray, masks_i[i])
ax[0][i].imshow(spectrum_img1, 'gray')
ax[0][i].set_xticks([])
ax[0][i].set_yticks([])
ax[1][i].imshow(img_ifft, 'gray')
ax[1][i].set_xticks([])
ax[1][i].set_yticks([])
ガウス分布型のハイパスフィルタ
fig, ax = plt.subplots(2, 3, figsize=(12, 8))
for i in range(3):
spectrum_img1, img_ifft, spectrum_img2 = masked_fft(img_gray, 1-g_masks[i])
ax[0][i].imshow(spectrum_img1/spectrum_img1.max(), 'gray', vmin=0, vmax=1)
ax[0][i].set_xticks([])
ax[0][i].set_yticks([])
ax[1][i].imshow(img_ifft, 'gray')
ax[1][i].set_xticks([])
ax[1][i].set_yticks([])
バンドパスフィルタ
中間的な周波数の範囲を残すようなものをバンドパスフィルタと呼びます。
マスク画像の作成をします。
masks_1 = np.array([np.array([make_mask(img_gray, r), make_mask(img_gray, r)]).transpose(1, 2, 0) for r in [30, 50, 100]])
masks_2 = np.array([np.array([make_mask(img_gray, r), make_mask(img_gray, r)]).transpose(1, 2, 0) for r in [50, 100, 180]])
mask_diff = masks_2 - masks_1
fig, ax = plt.subplots(2, 3, figsize=(12, 8))
for i in range(3):
spectrum_img1, img_ifft, spectrum_img2 = masked_fft(img_gray, mask_diff[i])
ax[0][i].imshow(spectrum_img1, 'gray')
ax[0][i].set_xticks([])
ax[0][i].set_yticks([])
ax[1][i].imshow(img_ifft, 'gray')
ax[1][i].set_xticks([])
ax[1][i].set_yticks([])
高域強調フィルタ
高域強調フィルタ
ハイパスフィルタは直流成分を含む低周波数成分を除去してしまいます。
低周波数成分を保ちつつ、高周波数成分を強調するフィルタを高域強調フィルタと呼びます。
以下のように、ハイパスフィルタを用いて求めることができます。
H_{h-emph}(u,v)=1+kH_{high}(u,v)
高域強調フィルタでは、エッジが強調され鮮鋭化された画像が得られます。
K = [1, 2, 3]
fig, ax = plt.subplots(2, 3, figsize=(12, 8))
for i in range(3):
spectrum_img1, img_ifft, spectrum_img2 = masked_fft(img_gray, K[i]+1-K[i]*g_masks[1])
ax[0][i].imshow(spectrum_img1, 'gray')
ax[0][i].set_xticks([])
ax[0][i].set_yticks([])
ax[1][i].imshow(img_ifft, 'gray')
ax[1][i].set_xticks([])
ax[1][i].set_yticks([])
空間フィルタリングによる鮮鋭化との関係
アンシャープマスキングと等価の処理となっている。
次回
幾何学的変換
参考