【Python】VGG16 を用いた Multilayer Stacked Covariance Pooling (MSCP)

ImageNet で事前学習した学習済みモデル (VGG16 など) を特徴量抽出器として，土地利用の画像分類 (UC Merced Land Use Dataset[^2] など) に使用する論文 ``Remote sensing scene classification using multilayer stacked covariance pooling[^1]'' の紹介と Python での実装．

# 実行環境

* Python 3.7.2
* Keras 2.4.3

# MSCP の概要

Multilayer Stacked Covariance Pooling (MSCP) のフレームワークを用いた画像分類は以下の手順で行われる．

1. 学習済みモデルのいくつかの隠れ層から多層の特徴量マップを抽出する．
1. 得られた特徴量マップを一つの多層特徴量マップにまとめ (stack)，まとめた特徴量マップから分散共分散行列を計算する．
1. 計算された分散共分散行列を特徴量として，線形カーネルを用いたサポートベクターマシン (SVM) を学習する．

今回は，学習済みモデルとして VGG16 を使用し，UC Merced Land Use Dataset[^2] を用いて画像分類を行う．  
データセットのディレクトリ構造は

```
UCMerced_LandUse
└── Images
    ├── agricultural
    |   ├── agricultural00.tif
    |   ├── ...
    |   └── agricultural99.tif
    ├── airplane
    ├── ...
```

のようになっている．

# 特徴量マップの抽出

Keras で VGG16 を使用するにあたり，以下のサイトを参考にした．

* [TensorFlow, KerasでVGG16などの学習済みモデルを利用 | note.nkmk.me](https://note.nkmk.me/python-tensorflow-keras-applications-pretrained-models/)
* [畳み込みニューラルネットワーク（CNN）の中間層の出力を可視化 | 人工知能（AI） | 有限会社はじめ研究所](https://hajimerobot.co.jp/ai/cnn_vgg16_layer/)

## VGG16 の層構造

Keras では，`keras.applications.vgg16.VGG16` で VGG16 のモデルを使用することができる．今回は，全結合層は使用しないので`include_top=False`，ImageNet で事前学習した重みを使用するため `weights='imagenet'` としている．

```python
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 の出力をアウトプットとした新しいモデルを定義する．

```python
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 のモデルに適した前処理を施すことができる．

```python
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()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/394475/84388382-810c-9427-b73e-80f4fbee02a3.png)

# 分散共分散行列の計算

## 特徴量マップのスタック

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$ としている．

```python
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 classification[^3]'' で Region Covariance Descriptor として提案された．この論文中では，画像の RGB 値や輝度値の勾配などをスタックし，分散共分散行列を計算している．概要は以下の記事にまとめている．  

* [【Python】Region Covariance: 分散共分散行列とコンピュータビジョン - Qiita](https://qiita.com/Mrrmm252/items/8e61da17b11011ddb661)

本論文で使用される特徴量は，Region Covariance で使われていた輝度値の勾配などの畳み込み演算を，畳み込みニューラルネットの隠れ層の出力で置き換えたものといえる．  

特徴量マップ $M \in \mathbb{R}^{s \times s \times D}$ から分散共分散行列は以下の手順で計算される．

1. 特徴量 $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)$ にベクトル化する．
1. 次式で分散共分散行列 $C \in \mathbb{R}^{D \times D}$ を計算する．

    ```math
    \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}
    ```
1. $C \leftarrow C + \epsilon I \ (\epsilon > 0)$ として，$C \in \mathbb{R}^{D \times D}$ の正定値性を保証する．

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/394475/26591808-35bf-110b-81cd-bed051484220.png)

Python では，`numpy.cov` で分散共分散行列を計算できる．

```python
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】行列指数関数・行列対数関数 - Qiita](https://qiita.com/Mrrmm252/items/a50a9b352e5064e40cc1)

Python では，エルミート行列用の固有値分解 `np.linalg.eigh` を用いて，行列の対数を実装できる．

```python
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, ...) に数字を割り当てる．

```python
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 を適用する準備ができたので，実際の流れを確認する．  

まず，ファイル名からラベルを，画像から分散共分散行列を抽出する．

```python
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)
```

次に，データを学習用と検証用に分割する．

```python
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 で学習・予測を行う．

```python
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 を確認する．

```python
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()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/394475/87bca528-cc57-b5a0-ec72-983004f13f64.png)


[^1]: [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.](https://www.umbc.edu/rssipl/people/aplaza/Papers/Journals/2018.TGRS.Multilayer.pdf)
[^2]: [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.](http://weegee.vision.ucmerced.edu/datasets/landuse.html)
[^3]: [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.](https://link.springer.com/content/pdf/10.1007/11744047_45.pdf)
[^4]: [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.](https://epubs.siam.org/doi/pdf/10.1137/050637996?casa_token=0qbcdZ1SRwwAAAAA:qxlEfIY0Krq1Oeef04H247PskamYd7HL-tnkxIdt8St2MWcrw0GjjpVTTQLoxs6aJExLW0V5eQeMYQ)