2. 画像フィルタ

本節では、OpenCVを用いた画像フィルタリングの方法について紹介する。

2.1. 単純なフィルタの実装

OpenCVを使った画像フィルタを行う前に、Pythonのfor文を使って画像フィルタを実装してみよう。

以下では、5✕5のGaussianフィルタを実装する。

# Gaussianカーネル
K = np.array(
    [
        [1, 4, 6, 4, 1],
        [4, 16, 24, 16, 4],
        [6, 24, 36, 24, 6],
        [4, 16, 24, 16, 4],
        [1, 4, 6, 4, 1],
    ],
    dtype=np.float32,
)
K /= K.sum()

次に適当な画像を読み込んで、フィルタリングを行う。この際、どのくらい時間がかかるのかを計測しておく。

height, width, _ = img.shape
print(f'Size: {width}x{height}')

start_time = time.time()
result = np.zeros_like(img)
for y in range(height):
    for x in range(width):
        color = np.zeros(3, dtype=np.float32)
        for dy in range(-2, 2 + 1):
            for dx in range(-2, 2 + 1):
                nx = np.clip(y + dy, 0, height - 1)
                ny = np.clip(x + dx, 0, width - 1)
                color += img[nx, ny] * K[dy + 2, dx + 2]

        result[y, x] = color

elapsed_time = time.time() - start_time
print(f'Elapsed time: {elapsed_time:.2f} sec')

Size: 512x512

Elapsed time: 63.46 sec

計算が完了したら、入力画像と出力画像を並べて表示してみよう。

fig = plt.figure()

ax1 = fig.add_subplot(1, 2, 1)
ax1.imshow(img)
ax1.set_title('Input')
ax1.axis('off')

ax2 = fig.add_subplot(1, 2, 2)
ax2.imshow(result)
ax2.set_title('Blurred')
ax2.axis('off')

fig.tight_layout()
plt.show()

残念なことだが、Pythonのfor文を用いたフィルタリングは非常に遅いことが分かる。そのため、OpenCVを始めとする画像処理用のライブラリはC++などの処理が高速な言語で実装されていることが多い。

実際に、OpenCVを使ったGaussianフィルタの例を見てみよう。フィルタに用いる関数は cv2.GaussianBlur() である。この関数は、入力画像、カーネルサイズ、Gauss関数の標準偏差 (= \(\sigma\)) を引数に取る。最後の \(\sigma\) に対して0が指定された場合には、カーネルサイズから自動的に \(\sigma\) が計算される。

start_time = time.time()
result = cv2.GaussianBlur(img, (5, 5), 0)
elapsed_time = time.time() - start_time
print(f'Elapsed time (OpenCV): {elapsed_time:.6f} sec')

Elapsed time (OpenCV): 0.001451 sec

fig = plt.figure()

ax1 = fig.add_subplot(1, 2, 1)
ax1.imshow(img)
ax1.set_title('Input')
ax1.axis('off')

ax2 = fig.add_subplot(1, 2, 2)
ax2.imshow(result)
ax2.set_title('Blurred')
ax2.axis('off')

fig.tight_layout()
plt.show()

このようにOpenCVに実装されたGaussianフィルタの処理は、Pythonで単純に実装したものと比べて非常に高速であることが分かる。

言うまでもなく、様々なフィルタがどのような原理で動作しているかを理解することは非常に重要である 。一方で、実際の応用においては、OpenCVなどのライブラリを利用して必要な処理を行い、他のより重要な処理に時間を割くことも重要だろう。

2.2. OpenCVにおける基本的なフィルタ

OpenCV上には、前述のGaussianフィルタ以外にも、様々な基本フィルタが実装されている。 画像フィルタ で紹介したフィルタの多くはOpenCV上で利用可能である。以下では、OpenCVにおける基本的なフィルタを紹介する。

img = skimage.data.coins()
img = img[:, :, None]

2.2.1. Sobelフィルタ

edge_x = cv2.Sobel(img, cv2.CV_8U, 0, 1)
edge_y = cv2.Sobel(img, cv2.CV_8U, 1, 0)

fig = plt.figure()

ax1 = fig.add_subplot(131)
ax1.imshow(img, cmap='gray')
ax1.set_title('Input')
ax1.axis('off')

ax2 = fig.add_subplot(132)
ax2.imshow(edge_x, cmap='gray')
ax2.set_title('Edge X')
ax2.axis('off')

ax3 = fig.add_subplot(133)
ax3.imshow(edge_y, cmap='gray')
ax3.set_title('Edge Y')
ax3.axis('off')

fig.tight_layout()
plt.show()

2.2.2. Cannyエッジ検出

edge = cv2.Canny(img, 100, 200)

fig = plt.figure()

ax1 = fig.add_subplot(121)
ax1.imshow(img, cmap='gray')
ax1.axis('off')

ax2 = fig.add_subplot(122)
ax2.imshow(edge, cmap='gray')
ax2.axis('off')

fig.tight_layout()
plt.show()

2.2.3. バイラテラル・フィルタ

img = skimage.data.coffee()
img = (img / 255.0).astype(np.float32)

res = cv2.bilateralFilter(img, -1, sigmaColor=0.3, sigmaSpace=11)

fig = plt.figure()

ax1 = fig.add_subplot(121)
ax1.imshow(img, cmap='gray')
ax1.set_title('Input')
ax1.axis('off')

ax2 = fig.add_subplot(122)
ax2.imshow(res, cmap='gray')
ax2.set_title('Bilateral Filtered')
ax2.axis('off')

fig.tight_layout()
plt.show()

2.2.4. 任意の線形フィルタ

線形フィルタをOpenCVを用いて実装する際には、カーネル関数を何らかの離散的な二次元テーブルとして用意したうえで cv2.filter2D を用いる。以下に3×3のGaussianフィルタの実装例を示す。

なお、以下のコードで cv2.filter2D の第2引数に与えている -1 は出力画像の型を表す変数で、入力画像を同じ型で出力する場合には -1 を指定すれば良い (cv2.Sobel の場合と同様)。

# フィルタカーネルの用意 (エンボス)
K = np.array(
    [
        [-2, -1, 0],
        [-1, 1, 1],
        [0, 1, 2],
    ],
    dtype=np.float32,
)

# フィルタの適用
res = cv2.filter2D(img, -1, K)
res = np.clip(res, 0, 1)

fig = plt.figure()

ax1 = fig.add_subplot(121)
ax1.imshow(img)
ax1.set_title('入力画像')
ax1.axis('off')

ax2 = fig.add_subplot(122)
ax2.imshow(res)
ax2.set_title('フィルタ結果')
ax2.axis('off')

plt.tight_layout()
plt.show()

2.2.5. 画像品質の評価

講義資料で解説したPSNRとSSIMを画像品質の評価に使うには、scikit-imageに実装された関数を使うと良い。

いずれも、 skimage.metrics モジュールに実装されており、PSNRは peak_signal_noise_ratio 関数、SSIMは structural_similarity 関数を用いる。両者とも、入力画像と出力画像を引数に取る。

使用上の注意点として、これらの画像を誤解なく使うのであれば、画像の値を8bit符号なし整数で表しておくのが良い。

また、SSIMはグレースケール画像に対して計算されるのが通例であり、カラー画像に対して用いる場合には、RGBの各チャンネルに対してSSIMを計算し、その平均を取る方法が一般的である。上記のscikit-imageの structural_similarity 関数は channel_axis 引数を適切に指定することで、カラー画像に対してもSSIMを計算できるようになっている。

import cv2
import matplotlib.pyplot as plt

import skimage
from skimage.metrics import peak_signal_noise_ratio as psnr
from skimage.metrics import structural_similarity as ssim

img = skimage.data.astronaut()
res1 = cv2.GaussianBlur(img, (5, 5), 0)
res2 = cv2.GaussianBlur(img, (9, 9), 0)

psnr1 = psnr(img, res1)
ssim1 = ssim(img, res1, channel_axis=2)

psnr2 = psnr(img, res2)
ssim2 = ssim(img, res2, channel_axis=2)

fig = plt.figure(figsize=(10, 4))

ax1 = fig.add_subplot(131)
ax1.imshow(img)
ax1.set_title('Input')
ax1.axis('off')
ax2 = fig.add_subplot(132)
ax2.imshow(res1)
ax2.set_title(f'GaussianBlur (5x5)\nPSNR: {psnr1:.2f} dB\nSSIM: {ssim1:.4f}')
ax2.axis('off')
ax3 = fig.add_subplot(133)
ax3.imshow(res2)
ax3.set_title(f'GaussianBlur (9x9)\nPSNR: {psnr2:.2f} dB\nSSIM: {ssim2:.4f}')
ax3.axis('off')

fig.tight_layout()
plt.show()

2.3. Fourier変換を用いたフィルタリング

2.3.1. NumPyによるFourier変換

NumPyを用いて二次元Fourier変換を行うには np.fft.fft2 を用いれば良い。この関数は axes=... に指定された2つの軸に沿ってFourier変換を行う。画像の場合には、第0軸が画像の行、第1軸が画像の列に対応するので、axes=(0, 1) と指定すれば良い。

なお、逆Fourier変換は np.fft.ifft2 を用いる。以下に、NumPyを用いたFourier変換の実装例を示す。

import numpy as np
import skimage.data

img = skimage.data.astronaut()  # 入力画像
img = (img / 255.0).astype(np.float32)  # float32に変換
freq = np.fft.fft2(img, axes=(0, 1))  # 2次元FFT

なお、Fourier変換の結果は複素数であるため、そのままでは表示ができない。そこで、各周波数成分の大きさをノルムとして求めて、その強度を表示してみる。このような各周波数における信号強度をパワースペクトルと呼ぶ。

パワースペクトルを画像として表示する際には、その対数をとって表示することが多い。パワースペクトルの表示例を以下に示す。

import matplotlib.pyplot as plt
import japanize_matplotlib

japanize_matplotlib.japanize()

# パワースペクトルの計算 (振幅の二乗)
ps = np.sum(np.abs(freq) ** 2, axis=2)

# パワースペクトル画像の描画
fig, ax = plt.subplots()
ax.imshow(np.log(ps + 1.0e-12), cmap='gray')  # 対数を取って表示
ax.set_title('パワースペクトル画像 (コーナーが原点)')
plt.show()

この図から分かるように、np.fft.fft2 の結果は画像のコーナーが原点になるような座標系を使っている。この座標系を画像の中心が原点になるような座標系に直すには np.fft.fftshift を用いる。

# 原点を中心にする
freq = np.fft.fftshift(freq, axes=(0, 1))

# パワースペクトルの計算 (振幅の二乗)
ps = np.sum(np.abs(freq) ** 2, axis=2)

# パワースペクトル画像の描画
fig, ax = plt.subplots()
ax.imshow(np.log(ps + 1.0e-12), cmap='gray')  # 対数を取って表示
ax.set_title('パワースペクトル画像 (中心が原点)')
plt.show()

ここから逆Fourier変換で元の画像を復元する際には、まず np.fft.ifftshift を用いて画像のコーナーが原点の座標系に戻してから、np.fft.ifft2 を用いて信号を画像空間に戻す。以下にその実装例を示す。

freq = np.fft.ifftshift(freq, axes=(0, 1))  # 原点を画像のコーナーに戻す
repro = np.fft.ifft2(freq, axes=(0, 1))  # 逆FFT
repro = np.real(repro)  # 実部だけを取り出す (虚部はほとんどゼロ)
repro = np.clip(repro, 0.0, 1.0)  # 0-1にクリッピング

# 元画像と復元画像の表示
fig = plt.figure()

ax = fig.add_subplot(121)
ax.imshow(img, vmin=0.0, vmax=1.0)
ax.set_title('元画像')
ax.axis('off')
ax = fig.add_subplot(122)
ax.imshow(repro, vmin=0.0, vmax=1.0)
ax.set_title('Fourier変換・逆変換で復元')
ax.axis('off')

plt.tight_layout()
plt.show()

2.3.2. 周波数空間でのフィルタ処理

GaussianフィルタやSobelフィルタなど、画像空間でカーネル関数が定義されるフィルタを周波数空間で適用する場合、画像とカーネル関数の両方をFourier変換してから、周波数空間で要素積を取り、その後で逆Fourier変換を行えば良いと述べた。

しかし、Fourier変換の結果得られる周波数空間の画像サイズは元の画像サイズと同じであるため、周波数空間でカーネル関数のFourier変換と要素積を取るためには、 カーネル関数を元の画像と同じ座標系の上で事前計算しておく必要がある 。

このようなカーネル関数の事前計算においては、画像の中心が原点となるような座標系を用い、その結果をコーナーが原点になるような座標系に修正してから、Fourier変換で周波数領域に変換する。

この方法の実装例を以下に示す。

import numpy as np
import skimage.data

img = skimage.data.astronaut()  # 入力画像
img = (img / 255.0).astype(np.float32)  # 0-1の範囲の実数に変換
H, W, _ = img.shape

# (H, W) の二次元グリッドを作成
ys, xs = np.mgrid[0:H, 0:W]

# 原点が中心に成るようにシフト
ys = ys - H // 2
xs = xs - W // 2

# Gaussianカーネルの計算
sigma = 2.0
kern = np.exp(-(xs**2 + ys**2) / (2 * sigma**2))
kern = kern / np.sum(kern)  # 正規化

# Fourier変換
K = np.fft.ifftshift(kern, axes=(0, 1))  # 原点を画像のコーナーに移動
K = np.fft.fft2(K, axes=(0, 1))  # FFT

# 画像もFourier変換
freq = np.fft.fft2(img, axes=(0, 1))  # FFT

# 周波数空間で要素積を取る
freq = freq * K[..., None]

# 逆FFT
res = np.fft.ifft2(freq, axes=(0, 1))  # 逆FFT
res = np.real(res)  # 実部だけを取り出す
res = np.clip(res, 0.0, 1.0)  # 0-1にクリッピング

# 元画像とフィルタ結果の表示
fig = plt.figure()
ax = fig.add_subplot(121)
ax.imshow(img, vmin=0.0, vmax=1.0)
ax.set_title('元画像')
ax.axis('off')
ax = fig.add_subplot(122)
ax.imshow(res, vmin=0.0, vmax=1.0)
ax.set_title('Gaussianフィルタ結果')
ax.axis('off')
plt.tight_layout()
plt.show()

上記の実装において、畳み込みカーネル K を変更することで、Gaussianフィルタ以外にもWienerフィルタなどの様々なフィルタが実装できる。

ただし、周波数領域でフィルタ処理を行う場合には画像のFourier変換と、カーネル関数のFourier変換の結果において、どこを原点とした座標系を用いているかには常に留意が必要である。

2.3.3. 画像の端での回り込み

上記のGaussianフィルタの結果をよく見てみると、画像の両端で不自然な模様が現れていることが分かる (特に右上の部分が分かりやすい)。

これは、Fourier変換を用いたフィルタリングにおいては、画像が上下左右にタイル状に繰り返されているように扱われるためで、画像の端では、同じ画像の反対側の端の情報が回り込んできてしまう。

これを防ぐには、画像の端の画素の値を使って、画像をパッディングしておくのが良い。このようなパッディング処理にはNumPyの np.pad が使用できる。例えば、端に10ピクセルのパッディングを行うためには、次のようにすればよい。

pad_img = np.pad(img, ((10, 10), (10, 10), (0, 0)), mode='edge')

この処理を行うことで、畳み込むカーネル関数のサイズにはよるものの、画像の回り込みの影響を減らすことができる。

2.4. 補足: Pythonでの画像フィルタの実装

Pythonで2次元のフィルタを実装する場合、素直にfor文を二重にしてしまうと、計算に時間がかかりすぎてしまう。これを避けるには、C++で実装すれば良い NumPyの機能を上手く使って実装する必要がある。

特に非線形フィルタを実装する場合には、単なる畳み込み演算にはならず、OpenCVに実装されている二次元畳み込みの関数 (cv2.filter2D) が利用できないので、何らかの工夫が必要になる。

2.4.1. カーネルとの重み付き和の計算をNumPyで行う

カーネル関数との畳み込みは各画素に対する二重ループ (x, yの各座標に対して1次元のループ)が必要で、さらに二次元のカーネル関数との画素の重み付き和の計算に二重ループが必要になる。しがって、全体としては四重のループが必要になり、for文の繰り返し計算に特に時間を要するPythonでは、比較的小さな画像であっても、かなりの時間を要する。

そこで、まず考えられる工夫はバイラテラル・フィルタのフィルタカーネルをNumPyの二次元配列に対する操作によって効率的に計算する、という方法である。

バイラテラル・フィルタのカーネル関数は色の差に関するGauss関数と空間的な距離に関するGauss関数の積であったので、次のような形で画像の一部を切り出すことでNumPyの配列操作により求めることが可能である。

# フィルタ範囲の取り出し
region = img_pad[y : y + ksize, x : x + ksize, :]
center = img_pad[y + pad, x + pad, :]
# 色カーネルの計算
d = region - center
d2 = np.sum(d * d, axis=-1)
L = np.exp(-0.5 * d2 / (sigma_r * sigma_r))
# カーネルの合成
W = K * L
W /= np.sum(W)
# 出力画素値の計算
res[y, x] = np.sum(region * W[..., None], axis=(0, 1))

少々、長文のコードになるが、以下に、この考え方を用いたバイラテラル・フィルタの実装例を示す。

Show code cell content

Hide code cell content

import time

import numpy as np
import skimage.data
import matplotlib.pyplot as plt
import japanize_matplotlib

# 画像の用意 (単精度浮動小数に変換しておく)
img = skimage.data.chelsea()
img = (img / 255.0).astype(np.float32)
height, width, _ = img.shape

# 時間計測開始
start_time = time.perf_counter()

# フィルタのパラメータ
ksize = 9
sigma_s = 5.0
sigma_r = 0.1

# 位置カーネルの事前計算
kx = np.arange(ksize) - ksize // 2
ky = np.arange(ksize) - ksize // 2
ky, kx = np.meshgrid(ky, kx, indexing='ij')
K = np.exp(-0.5 * (kx * kx + ky * ky) / (sigma_s * sigma_s))

# 画像をあらかじめパディング
pad = ksize // 2
img_pad = np.pad(img, ((pad, pad), (pad, pad), (0, 0)), mode='edge')

# フィルタ処理
res = np.zeros_like(img)
for y in range(height):
    for x in range(width):
        # フィルタ範囲の取り出し
        region = img_pad[y : y + ksize, x : x + ksize, :]
        center = img_pad[y + pad, x + pad, :]
        # 色カーネルの計算
        d = region - center
        d2 = np.sum(d * d, axis=-1)
        L = np.exp(-0.5 * d2 / (sigma_r * sigma_r))
        # カーネルの合成
        W = K * L
        W /= np.sum(W)
        # 出力画素値の計算
        res[y, x] = np.sum(region * W[..., None], axis=(0, 1))

# 時間計測終了
end_time = time.perf_counter()
elapsed_time = end_time - start_time

# 画像の表示
fig = plt.figure()
gs = fig.add_gridspec(1, 2)
ax = plt.subplot(gs[0])
ax.imshow(img)
ax.set_title('元画像')
ax.axis('off')
ax = plt.subplot(gs[1])
ax.imshow(res)
ax.set_title(f'バイラテラル・フィルタ ({elapsed_time:.2f}秒)')
ax.axis('off')
fig.tight_layout()
plt.show()

2.4.2. スライドウィンドウ分割の利用

上記の方法でも、四重ループのfor文を用いるよりはかなり計算が早くなるが、それでも450×300の画像に対して1秒程度要しており、決して高速とは言えない。

ここまでのコードでは画像の各画素の周囲の領域を取り出すためにNumPyのスライス (img_pad[y:y+ksize, x:x+ksize]のような書き方) を用いたが、この画像の一部の切り出しもfor文を用いる代わりにNumPyの関数 (C,C++等の言語で実装されており早い) を用いれば、より高速に計算できる可能性がある。

実際、NumPyには、numpy.lib.stride_tricks というモジュール内に slinding_window_view という関数が用意されていて、\((H, W, 3)\) という形の画像 img_pad に対して、 sliding_window_view(img_pad, (ksize, ksize, 3)) のように呼び出すと、画像が (H, W, ksize, ksize, 3) という形の5次元の配列に変換される。

この結果を用いれば、前述のカーネル関数の計算をすべての画素に対してNumPyの配列処理によって計算できるようになるので、さらなる高層化が期待できる。

再び、長文のコードになるが、以下に sliding_window_view を用いたバイラテラル・フィルタの実装例を示す。

Show code cell content

Hide code cell content

import time

import numpy as np
import skimage.data
import matplotlib.pyplot as plt
import japanize_matplotlib
from numpy.lib.stride_tricks import sliding_window_view

# 画像の用意 (単精度浮動小数に変換しておく)
img = skimage.data.chelsea()
img = (img / 255.0).astype(np.float32)
height, width, _ = img.shape

# 時間計測開始
start_time = time.perf_counter()

# フィルタのパラメータ
ksize = 9
sigma_s = 5.0
sigma_r = 0.1

# 位置カーネルの事前計算
kx = np.arange(ksize) - ksize // 2
ky = np.arange(ksize) - ksize // 2
ky, kx = np.meshgrid(ky, kx, indexing='ij')
K = np.exp(-0.5 * (kx * kx + ky * ky) / (sigma_s * sigma_s))

# 画像をあらかじめパディング
pad = ksize // 2
img_pad = np.pad(img, ((pad, pad), (pad, pad), (0, 0)), mode='edge')

# 画素ごとのフィルタに必要な画像範囲をパッチとして取り出す
img_view = sliding_window_view(img_pad, (ksize, ksize, 3))
img_view = img_view.reshape(img.shape[0], img.shape[1], ksize, ksize, 3)
center = img_view[:, :, pad : pad + 1, pad : pad + 1, :]

# 色カーネルの計算
diff = img_view - center
d2 = np.sum(diff * diff, axis=-1)
L = np.exp(-0.5 * d2 / (sigma_r * sigma_r))

# カーネルの合成
W = K * L
W = W / np.sum(W, axis=(-1, -2), keepdims=True)
W = W.reshape(img.shape[0], img.shape[1], ksize, ksize, 1)

# 各パッチを一気に畳み込む
res = np.sum(img_view * W, axis=(-2, -3))

# 時間計測終了
end_time = time.perf_counter()
elapsed_time = end_time - start_time

# 画像の表示
fig = plt.figure()
gs = fig.add_gridspec(1, 2)
ax = plt.subplot(gs[0])
ax.imshow(img)
ax.set_title('元画像')
ax.axis('off')
ax = plt.subplot(gs[1])
ax.imshow(res)
ax.set_title(f'バイラテラル・フィルタ ({elapsed_time:.2f}秒)')
ax.axis('off')
fig.tight_layout()
plt.show()

こちらの実装は、この前に示した実装と比べて2-3倍程度早くなっており、450×300の画像に対して0.5秒と、まだまだ満足行く測度ではないもののPythonに閉じた実装としては、かなり高速に計算できるようになっている。

なお、実際にバイラテラル・フィルタが使えれば良いだけであれば、OpenCVに実装されている cv2.bilateralFilter を用いるのが高速かつ簡単だろう。