ImageNet で事前学習した学習済みモデル (VGG16 など) を特徴量抽出器として,土地利用の画像分類 (UC Merced Land Use Dataset1 など) に使用する論文 ``Remote sensing scene classification using multilayer stacked covariance pooling2'' の紹介と Python での実装.
実行環境
- Python 3.7.2
- Keras 2.4.3
MSCP の概要
Multilayer Stacked Covariance Pooling (MSCP) のフレームワークを用いた画像分類は以下の手順で行われる.
- 学習済みモデルのいくつかの隠れ層から多層の特徴量マップを抽出する.
- 得られた特徴量マップを一つの多層特徴量マップにまとめ (stack),まとめた特徴量マップから分散共分散行列を計算する.
- 計算された分散共分散行列を特徴量として,線形カーネルを用いたサポートベクターマシン (SVM) を学習する.
今回は,学習済みモデルとして VGG16 を使用し,UC Merced Land Use Dataset1 を用いて画像分類を行う.
データセットのディレクトリ構造は
UCMerced_LandUse
└── Images
├── agricultural
| ├── agricultural00.tif
| ├── ...
| └── agricultural99.tif
├── airplane
├── ...
のようになっている.
特徴量マップの抽出
Keras で VGG16 を使用するにあたり,以下のサイトを参考にした.
- TensorFlow, KerasでVGG16などの学習済みモデルを利用 | note.nkmk.me
- 畳み込みニューラルネットワーク(CNN)の中間層の出力を可視化 | 人工知能(AI) | 有限会社はじめ研究所
VGG16 の層構造
Keras では,keras.applications.vgg16.VGG16
で VGG16 のモデルを使用することができる.今回は,全結合層は使用しないのでinclude_top=False
,ImageNet で事前学習した重みを使用するため weights='imagenet'
としている.
from keras.applications.vgg16 import VGG16
vgg16 = VGG16(include_top=False, weights='imagenet')
vgg16.summary()
# Model: "vgg16"
# _________________________________________________________________
# Layer (type) Output Shape Param #
# =================================================================
# input_2 (InputLayer) [(None, None, None, 3)] 0
# _________________________________________________________________
# block1_conv1 (Conv2D) (None, None, None, 64) 1792
# _________________________________________________________________
# block1_conv2 (Conv2D) (None, None, None, 64) 36928
# _________________________________________________________________
# block1_pool (MaxPooling2D) (None, None, None, 64) 0
# _________________________________________________________________
# block2_conv1 (Conv2D) (None, None, None, 128) 73856
# _________________________________________________________________
# block2_conv2 (Conv2D) (None, None, None, 128) 147584
# _________________________________________________________________
# block2_pool (MaxPooling2D) (None, None, None, 128) 0
# _________________________________________________________________
# block3_conv1 (Conv2D) (None, None, None, 256) 295168
# _________________________________________________________________
# block3_conv2 (Conv2D) (None, None, None, 256) 590080
# _________________________________________________________________
# block3_conv3 (Conv2D) (None, None, None, 256) 590080
# _________________________________________________________________
# block3_pool (MaxPooling2D) (None, None, None, 256) 0
# _________________________________________________________________
# block4_conv1 (Conv2D) (None, None, None, 512) 1180160
# _________________________________________________________________
# block4_conv2 (Conv2D) (None, None, None, 512) 2359808
# _________________________________________________________________
# block4_conv3 (Conv2D) (None, None, None, 512) 2359808
# _________________________________________________________________
# block4_pool (MaxPooling2D) (None, None, None, 512) 0
# _________________________________________________________________
# block5_conv1 (Conv2D) (None, None, None, 512) 2359808
# _________________________________________________________________
# block5_conv2 (Conv2D) (None, None, None, 512) 2359808
# _________________________________________________________________
# block5_conv3 (Conv2D) (None, None, None, 512) 2359808
# _________________________________________________________________
# block5_pool (MaxPooling2D) (None, None, None, 512) 0
# =================================================================
# Total params: 14,714,688
# Trainable params: 14,714,688
# Non-trainable params: 0
# _________________________________________________________________
VGG16 の隠れ層の取得
論文中では,block3_conv3, block4_conv3, block5_conv3 の出力を特徴量マップとして使用している.そこで,block3_conv3, block4_conv3, block5_conv3 の出力をアウトプットとした新しいモデルを定義する.
from keras.models import Model
outputs = [vgg16.layers[i].output for i in (9, 13, 17)]
model = Model(inputs=vgg16.input, outputs=outputs)
実際に隠れ層の出力を確認してみる.入力画像には,keras.applications.vgg16.preprocess_input
で VGG16 のモデルに適した前処理を施すことができる.
from typing import List
from keras.applications.vgg16 import preprocess_input
import matplotlib.pyplot as plt
import numpy as np
from PIL import Image
def extract_feature(model: Model,
img: np.ndarray) -> List[np.ndarray]:
# reshape (H, W, 3) -> (1, H, W, 3)
processed_img = np.expand_dims(img, axis=0)
processed_img = preprocess_input(processed_img)
preds = model.predict(processed_img)
# reshape (1, H, W, 3) -> (H, W, 3)
preds = [pred[0] for pred in preds]
return preds
file = './UCMerced_LandUse/Images/river/river00.tif'
img = np.array(Image.open(file))
features = extract_feature(model, img)
n_cols = len(features)+1
plt.figure(figsize=(16, 4))
plt.subplot(1, n_cols, 1)
plt.imshow(img)
plt.title(f'Input\nshape={img.shape}')
plt.axis('off')
for i, feature in enumerate(features, start=2):
plt.subplot(1, n_cols, i)
plt.imshow(feature[:, :, 0], cmap='jet')
plt.title(f'block{i+1}_conv3\nshape={feature.shape}')
plt.axis('off')
plt.show()
分散共分散行列の計算
特徴量マップのスタック
block3_conv3, block4_conv3, block5_conv3 の出力をそれぞれ,$M_{3,3}, M_{4,3}, M_{5,3}$ とすると,上で見たように,$M_{3,3} \in \mathbb{R}^{64 \times 64 \times 256}, M_{4,3} \in \mathbb{R}^{32 \times 32 \times 512}, M_{5,3} \in \mathbb{R}^{16 \times 16 \times 512}$ である.これらの特徴量マップをスタックするために,バイリニア補間を用いて $s \times s$ にリサイズする.
さらに,特徴量マップのチャネル数 (今回の場合は 256, 512, 512) の和が計算する分散共分散行列の次元数となるので,チャネルごとの平均を取ることにより,$d$ 次元まで次元削減する.
結果,$M_{3,3}, M_{4,3}, M_{5,3}$ はそれぞれ $\hat M_{3,3}, \hat M_{4,3}, \hat M_{5,3} \in \mathbb{R}^{s \times s \times d}$ へとダウンサンプルされ,これらをスタックすることにより,
$$M = \left[\hat M_{3,3};\hat M_{4,3};\hat M_{5,3}\right] \in \mathbb{R}^{s \times s \times D}, \quad D = 3d$$
が特徴量マップとして得られる.
ここまでを Python でまとめると以下の通り.今回は $s=16, d=32$ としている.
import cv2
def stack(features: List[np.ndarray],
depth: int = 32, size: int = 16) -> np.ndarray:
features = [cv2.resize(feature, (size, size)) for feature in features]
features = np.dstack(
[np.dstack([np.mean(splitted, axis=-1)
for splitted in np.array_split(feature, depth, axis=-1)])
for feature in features])
return features
depth = 32
size = 16
features = stack(features, depth=depth, size=size)
features.shape
# (16, 16, 96)
分散共分散行列の計算
分散共分散行列を画像の特徴量として使用するというアイデアは,2006 年に ``Region covariance: A fast descriptor for detection and classification3'' で Region Covariance Descriptor として提案された.この論文中では,画像の RGB 値や輝度値の勾配などをスタックし,分散共分散行列を計算している.概要は以下の記事にまとめている.
本論文で使用される特徴量は,Region Covariance で使われていた輝度値の勾配などの畳み込み演算を,畳み込みニューラルネットの隠れ層の出力で置き換えたものといえる.
特徴量マップ $M \in \mathbb{R}^{s \times s \times D}$ から分散共分散行列は以下の手順で計算される.
-
特徴量 $M \in \mathbb{R}^{s \times s \times D}$ を $[z_1, z_2, \ldots, z_N] \in \mathbb{R}^{N \times D}\ (N = s^2)$ にベクトル化する.
-
次式で分散共分散行列 $C \in \mathbb{R}^{D \times D}$ を計算する.
\begin{align} \mu &= \frac1N \sum_{i=1}^N z_i \in \mathbb{R}^D \\ C &= \frac{1}{N-1} \sum_{i=1}^N (z_i - \mu) (z_i - \mu)^\top \in \mathbb{R}^{D \times D} \end{align}
-
$C \leftarrow C + \epsilon I \ (\epsilon > 0)$ として,$C \in \mathbb{R}^{D \times D}$ の正定値性を保証する.
Python では,numpy.cov
で分散共分散行列を計算できる.
epsilon = 1e-8
cov = np.cov(features.reshape(-1, depth * 3), rowvar=False)
cov += np.identity(depth * 3) * epsilon
SVM の適用
ユークリッド空間への変換
計算された分散共分散行列 $C \in \mathbb{R}^{D \times D}$ は正定値対称行列であり,ユークリッド空間の元ではないため,SVM をそのまま適用することができない.そこで,行列の対数を作用させることにより,ユークリッド空間の元 (対称行列) とする4.
$$\hat C = \mathop{\mathrm{logm}}C = U \log(\Sigma) U^\top \in \mathbb{R}^{D \times D}$$
ここで,$C=U\Sigma U^\top$ は分散共分散行列 $C$ の固有値分解である.行列の対数については,以下にまとめている.
Python では,エルミート行列用の固有値分解 np.linalg.eigh
を用いて,行列の対数を実装できる.
def logm(x: np.ndarray) -> np.ndarray:
w, v = np.linalg.eigh(x)
return v @ (np.log(w) * v).T
cov = logm(cov)
ラベルの準備
クラス名 (agricultural, airplane, ...) に数字を割り当てる.
import os
img_path = './UCMerced_LandUse/Images/'
class_names = sorted(os.listdir(img_path))
class_names = [class_name for class_name in class_names
if os.path.isdir(os.path.join(img_path, class_name))]
class_name_dict = {class_name: i for i, class_name in enumerate(class_names)}
class_name_dict
# {'agricultural': 0,
# 'airplane': 1,
# ...
# 'tenniscourt': 20}
SVM の適用
SVM を適用する準備ができたので,実際の流れを確認する.
まず,ファイル名からラベルを,画像から分散共分散行列を抽出する.
import glob
from tqdm.notebook import tqdm
# Hyperparameters
depth = 32
size = 16
epsilon = 1e-8
# Prepare model
vgg16 = VGG16(include_top=False, weights='imagenet')
outputs = [vgg16.layers[i].output for i in (9, 13, 17)]
model = Model(inputs=vgg16.input, outputs=outputs)
filenames = sorted(glob.glob('./UCMerced_LandUse/Images/*/*.tif'))
X, y = [], []
for file in tqdm(filenames):
class_name = os.path.basename(os.path.dirname(file))
y.append(class_name_dict[class_name])
# Load an image
img = np.array(Image.open(file))
# Extract feature maps
features = extract_feature(model, img)
# Stack
features = stack(features, depth=depth, size=size)
# Covariance
cov = np.cov(features.reshape(-1, depth * 3), rowvar=False)
cov += np.identity(depth * 3) * epsilon
# Logarithm
cov = logm(cov)
X.append(cov.ravel())
X = np.vstack(X)
y = np.array(y)
次に,データを学習用と検証用に分割する.
from sklearn.model_selection import train_test_split
train_X, test_X, train_y, test_y = train_test_split(
X, y, train_size=0.8, random_state=0, stratify=y)
SVM で学習・予測を行う.
from sklearn.svm import SVC
svm = SVC(C=1, kernel='linear', random_state=0).fit(train_X, train_y)
pred_y = svm.predict(test_X)
Accuracy, Confusion Matrix を確認する.
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.metrics import accuracy_score, confusion_matrix
print(f'Accuracy: {accuracy_score(test_y, pred_y)}')
# Accuracy: 0.9523809523809523
plt.figure(figsize=(8, 6))
sns.heatmap(confusion_matrix(test_y, pred_y), annot=True, square=True,
xticklabels=class_names, yticklabels=class_names)
plt.xlabel('Predicted')
plt.ylabel('Actual')
plt.show()
-
Yi Yang and Shawn Newsam, "Bag-Of-Visual-Words and Spatial Extensions for Land-Use Classification," ACM SIGSPATIAL International Conference on Advances in Geographic Information Systems (ACM GIS), 2010. ↩ ↩2
-
He, Nanjun, et al. "Remote sensing scene classification using multilayer stacked covariance pooling." IEEE Transactions on Geoscience and Remote Sensing 56.12 (2018): 6899-6910. ↩
-
Tuzel, Oncel, Fatih Porikli, and Peter Meer. "Region covariance: A fast descriptor for detection and classification." European conference on computer vision. Springer, Berlin, Heidelberg, 2006. ↩
-
Arsigny, Vincent, et al. "Geometric means in a novel vector space structure on symmetric positive-definite matrices." SIAM journal on matrix analysis and applications 29.1 (2007): 328-347. ↩