国土地理院の空中写真データセットでセマンティック・セグメンテーションしてみた

# はじめに

この記事では，国土地理院が公開している **セマンティック・セグメンテーション** 用の **空中写真データセット** を利用して，空中写真から道路領域を検出することに挑戦します．国土地理院の報告では，評価指標であるF値（Dice係数）が **0.8以上** となったそうです．このF値（Dice係数）と同等の精度を得ることを目標にPyTorchで実装して確認します．

# 地物抽出用の空中写真データセット

[空中写真データセット（GSIデータセット）](https://gisstar.gsi.go.jp/gsi-dataset/index_ja.html)が国土地理院より公開されました． この空中写真（航空写真）は，平成30年度〜令和4年度にかけて実施された国土地理院の特別研究において作成された機械学習用のデータセットです．後述の セマンティック・セグメンテーションと呼ばれる手法の訓練データとして利用が可能であり，国土地理院技術資料として広く公開することになったそうです．現在は **道路** ， **水部** ， **防波堤** ， **太陽光発電設備** ， **タンク** ， **透過水制** の6種類が公開されており，近日中に **普通建物** ， **茶畑** の2種類が公開されると発表されています．

| 対象 | 着色 |
|:-:|:-:|
| 道路  | 道路を赤色に着色 #FF0000 |
| 水部  | 水部を青色に着色 #0000FF |
| 防波堤  | 防波堤を黄色に着色 #FFFF00 |
| 太陽光発電設備  | 太陽光発電設備を黄色に着色 #FFFF00 |
| タンク  |  タンクを赤色に着色 #FF0000 |
| 透過水制  | 透過水制を黄色に着色 #FFFF00|

ここでは，[道路の空中写真データセット](https://gisstar.gsi.go.jp/gsi-dataset/02/index.html)に着目し，機械学習ライブラリの[PyTorch](https://pytorch.org/)を利用して，セマンティック・セグメンテーションを実装し，その精度を確認することを目的とします．道路のデータセットは， **286x286ピクセル** と **572x572ピクセル** の2種類の空中写真が提供されています．ここでは，前者の286x286ピクセルの空中写真を利用します．サンプル数は$N=10000$です．空中写真は**地上画素寸法20cm**で撮影されており，これは1ピクセルの長さが20センチメートルであることを意味しています．このため，画像の1辺の長さは57.2メートルです（次の式で算出）．

$$
286(px) \times 20(cm)=5720(cm)=57.2(m)
$$

データセットに含まれる空中写真の例を示します（上から順に32.png，240.png，615.png）．左がオリジナルの画像（空中写真），右が道路のマスク画像です．この色付け（ラベリング）は人手によって行われたそうです．マスク画像に含まれる道路領域が赤で着色されていることが確認できますが，田んぼの畦道や駐車場など，道路との境界が曖昧なケースもあるようです．

[![Image from Gyazo](https://i.gyazo.com/d6bc6e731be84f675b9dd1f5595acd71.png)](https://gyazo.com/d6bc6e731be84f675b9dd1f5595acd71)
<!-- 32.png -->

[![Image from Gyazo](https://i.gyazo.com/4b98e8f7ec9132eb24b0cad79a3be167.png)](https://gyazo.com/4b98e8f7ec9132eb24b0cad79a3be167)
<!-- 240.png -->

[![Image from Gyazo](https://i.gyazo.com/4ab7b753a91e9fd554b5966f32fe4262.png)](https://gyazo.com/4ab7b753a91e9fd554b5966f32fe4262)
<!-- 615.png -->

国土地理院のサイトには，このデータセットを利用することで，評価指標である **F値（Dice係数）** が0.8以上になったと記載されています．この値にどこまで近付けるか挑戦してみます！

# セマンティック・セグメンテーション

セマンティック・セグメンテーション（Semantic Segmentation）は，画像の画素をラベリング（分類とも考えられる）する手法のことです．セマンティック・セグメンテーション用のデータセットとして **COCO** や **CityScapes** などが有名です． COCOは，人（person），バックパック（backpack），傘（umbrella）など80種類のオブジェクトを対象としたデータセットです．また，CityScapesは，道路（road），車（car），建物（building）など都市の風景を対象としたデータセットです．

- [COCO(Common Objects in Context)](https://cocodataset.org/#home)
- [CityScapes](https://www.cityscapes-dataset.com/)

COCOデータセットに含まれるサンプル画像を次に示します．最初の画像には，人物（person），犬（dog），自転車（bicycle）などのラベルが色付けされています．2番目の画像には，人物（person），鳥（bird），ハンドバッグ（handbag）などのラベルが色付けされています．本記事で取り上げる空中写真データセットも，これらと同種のデータセットとみなすことができます．

[![Image from Gyazo](https://i.gyazo.com/a9abb2f704e2abfadee93ad01a85ee57.png)](https://gyazo.com/a9abb2f704e2abfadee93ad01a85ee57)

[![Image from Gyazo](https://i.gyazo.com/5885037a0e53677c2fdb83fda6f63d9f.png)](https://gyazo.com/5885037a0e53677c2fdb83fda6f63d9f)

## 学習モデル

セマンティック・セグメンテーションには，畳み込みニューラルネットワーク（CNN: Convolutional Neural Network）を拡張した[SegNet](https://arxiv.org/abs/1511.00561)，[U-Net](https://arxiv.org/abs/1505.04597)などのモデルが用いられます．
標準的なCNNは次の図で示すような構造となっています．入力画像は，畳み込み層＆プーリング層でダウン・サンプリングされることで特徴抽出がなされ，最後は全結合層で分類結果を表す **ラベル（人，犬など）** が出力されます．

[![Image from Gyazo](https://i.gyazo.com/92d87041c3d0c7c4aa11a93ed4c9e44a.png)](https://gyazo.com/92d87041c3d0c7c4aa11a93ed4c9e44a)

これに対し，SegNetは次の図で示すような構造となっています．Encoder-Decoder型と呼ばれる構造であり，入力画像はエンコーダでダウン・サンプリングされることで特徴抽出されます．その後，デコーダでダウン・サンプリングと同じだけのアップ・サンプリングされることで，入力画像に対応する **マスク画像** を出力します．U-NetもEncoder-Decoder型のモデルであり，**スキップ接続** と呼ばれる方法で，ダウン・サンプリングで消失される情報を保持する工夫が加えられています．後述のセマンティック・セグメンテーションの実装では，U-Netを採用します．

[![Image from Gyazo](https://i.gyazo.com/a79255cfb2e9911e27ffe99710c78a1e.png)](https://gyazo.com/a79255cfb2e9911e27ffe99710c78a1e)

## 評価指標

学習モデルの評価には **Dice Score** と **IoU Score**を採用することにします．また，学習の損失関数には，1からDice Scoreを引いた **Dice Loss**を用います．下図の例で考えてみましょう．左から「正解のマスク画像（訓練データ）」，「予測されたマスク画像」，「正しく予測されたピクセル（正解と予測が重なっているピクセル）」を表しています．画像は$5 \times 5 = 25$のピクセルで構成されています． 

[![Image from Gyazo](https://i.gyazo.com/cc0bd148ffe646d1f0b1007bf3b87199.png)](https://gyazo.com/cc0bd148ffe646d1f0b1007bf3b87199)

この画像の中で下記の条件を満たすピクセルを数えることで， **TP(True Positive)** ，**TN(True Negative)** ，**FP(False Positive)** ，**FN(False Negative)** を算出します．これらの値から後述の **適合率（Precision）** と **再現率（Recall）** が求められます．

- TP: マスクと予測し，実際にマスクだったピクセル数
- TN: 背景と予測し，実際に背景だったピクセル数
- FP: マスクと予測し，実際は背景だったピクセル数
- FN: 背景と予測し，実際はマスクだったピクセル数

|  | Positive | Negative |
|:-:|:-:|:-:|
| **True**  | TP=4  | TN=15  |
| **Fase**  | FP=4 | FN=2  |

### 適合率（Precision）

適合率は，マスクと予測されたピクセル数（$TP+FP$）に占める，実際にマスクだったピクセル数（$TP$）の割合です．

$$
Precision = \frac{TP}{TP+FP}=\frac{4}{4+4} = 0.5
$$

### 再現率（Recall）

再現率は，実際のマスクのピクセル数（$TP+FN$）に占める，マスクと予測されたピクセル数（$TP$）の割合です．

$$
Recall = \frac{TP}{TP+FN}=\frac{4}{4+2} \simeq 0.67
$$

### Dice Score（Dice係数）

F値とも呼ばれます．上記の適合率と再現率の調和平均です．適合率と再現率のバランスを考慮した評価指標です．予測されたマスクが正解のマスクと完全に一致したときに$Dice=1$になります．

$$
Dice = \frac{2}{\frac{1}{Precision} + \frac{1}{Recall}} = \frac{TP}{TP+\frac{FP+FN}{2}}= \frac{4}{4+\frac{4+2}{2}} \simeq 0.57
$$

また，損失関数として用いる **Dice Loss** は次の式で算出されます． Dice Lossの最小化は，Dice Scoreの最大化と同じ結果となります．

$$
DiceLoss = 1 - Dice \simeq 0.43
$$

### IoU Score (Intersection over Union Score）

ジャッカード係数（Jaccard Score）とも呼ばれます．正解のマスクのピクセル集合を$A$，予測されたマスクのピクセル集合を$B$とします．このとき，和集合$A \cup B$に占める，共通集合$A \cap B$の割合です．Dice Scoreと同様に，予測されたマスクが正解のマスクと完全に一致したときに$IoU=1$になります．

$$
IoU = \frac{TP}{TP+FP+FN}=\frac{4}{4+4+2}=0.4
$$

# 実装

[PyTorch](https://pytorch.org/)を利用して，空中写真のデータセットを学習するための[U-Net](https://arxiv.org/abs/1505.04597)を実装します．実装には[Segmentation Models](https://github.com/qubvel/segmentation_models.pytorch)というライブラリを利用します．

```python
# Segmentaion Modelsのインストール
!pip install -U segmentation-models-pytorch
```

空中写真のデータセットに含まれる画像のリストを次の形式の[pandas](https://pandas.pydata.org/)のデータフレーム`dataset_df`で保持することにします．ここで，　**org** はオリジナル画像のパス， **val** はマスク画像のパスです．データセットに含まれるサンプル数は$N=10000$ですが，ここでは簡単化のため，1.pngから1000.pngまでの$N＝1000$のサンプルを対象とします．

| index |id|org|val|
|:----|:----|:----|:----|
|0|map-1|./dataset/gsi-1000/org286/1.png|./dataset/gsi-1000/val286/1.png|
|1|map-2|./dataset/gsi-1000/org286/2.png|./dataset/gsi-1000/val286/2.png|
|2|map-3|./dataset/gsi-1000/org286/3.png|./dataset/gsi-1000/val286/3.png|
|3|map-4|./dataset/gsi-1000/org286/4.png|./dataset/gsi-1000/val286/4.png|
|4|map-5|./dataset/gsi-1000/org286/5.png|./dataset/gsi-1000/val286/5.png|

```python
 # データセットのデータフレーム
dataset_df = pd.read_csv(dataset_file)
```

データセットに含まれるマスク画像は，道路が赤（#ff0000）でマスクされています．このままでは，PyTorchのマスクとして利用することができません．そこで，次の`generateMask`関数を定義し，道路は`1`，その他は`0`のバイナリ画像（[Pillow](https://pillow.readthedocs.io/en/stable/)のImageオブジェクト）に変換します．

```python
# マスク画像の生成
def generateMask(image):

  width, height = image.size
  image_mask = Image.new("1", (width, height)) # "1"はバイナリ画像

  for x in range(width):
    for y in range(height):
      r, g, b = image.getpixel((x, y))
        
      # 赤のとき
      if(r == 255 and g == 0 and b == 0):
        image_mask.putpixel((x, y), (1)) # ピクセルの値を1に設定
      else:
        image_mask.putpixel((x, y), (0)) # ピクセルの値を0に設定

  return image_mask
```

PyTorchでデータセットを読み込みためのクラス`GSIDataset`を定義します．`__getitem__`関数は指定されたインデックスのサンプルを返します．サンプルはオリジナル画像`image_org`とマスク画像`image_mask`のタプルで表されます．また，`__len__`関数はデータセットに含まれるサンプル数を返します．

```python
# データを読み込むためのクラス
class GSIDataset(Dataset):

  def __init__(self, dataset_df, transform):
    self.dataset_df = dataset_df
    self.transform = transform
    self.id_list = dataset_df["id"].to_list()
    self.org_list = dataset_df["org"].to_list()
    self.var_list = dataset_df["val"].to_list()

  # サンプルの取得（オリジナル画像とマスク画像）
  def __getitem__(self, index):
    with open(self.org_list[index], "rb") as f_org:
      image_org = Image.open(f_org)
      image_org = self.transform(image_org)
    
    with open(self.var_list[index], "rb") as f_var:
      image_var = Image.open(f_var)
      image_mask = generateMask(image_var)　# マスク画像をPyTorch用に変換
      image_mask = self.transform(image_mask)

    return (image_org, image_mask)

  # サンプル数の取得
  def __len__(self):
    length = len(self.dataset_df)
    return length
```

交差検証のためデータセットを訓練用と評価用で$9:1$に分割します．訓練用のサンプル数は900，評価用のサンプル数は100となります．`train_transform`と`eval_transform`は画像に対する前処理です．ここでは，画像サイズを$286 \times 286$から，$256 \times 256$に縮小します（32の倍数という制約があるため）．加えて，画像をランダムに上下左右に反転することでデータ拡張（Data Augumentation）します．

```python
# データセットを訓練用と評価用に分割（9:1）
train_df, eval_df = train_test_split(dataset_df, shuffle=True, test_size=0.1, random_state=random_state)

train_df = train_df.reset_index(drop=True) # 訓練用
eval_df = eval_df.reset_index(drop=True) # 評価用

train_dataset = GSIDataset(train_df, train_transform) # 訓練用
eval_dataset = GSIDataset(eval_df, eval_transform) # 評価用
```

ミニバッチ学習用のデータローダを作成します．バッチサイズは`16`とします．

```python
# 学習用データローダ
train_dataloader = DataLoader(train_dataset, batch_size=16, shuffle=True)

# 評価用データローダ
eval_dataloader = DataLoader(eval_dataset, batch_size=16, shuffle=False)
```

U-Netのモデルを定義します．エンコーダには[ResNet](https://arxiv.org/abs/1512.03385)を採用します．入力はRGB画像であるため`in_channels=3`，1クラスの分類（"道路"と"その他"）のため`classes=1`，出力画像の各ピクセルは0から1の範囲で表すため`activation="sigmoid"`を設定します．

```python
# U-Net
network = smp.Unet(
    encoder_name="resnet34",
    encoder_weights="imagenet",
    in_channels = 3,
    classes = 1,
    activation = "sigmoid",
)
```

損失関数に **DiceLoss** ，評価指標に **IoU** ，最適化手法に[Adam](https://arxiv.org/abs/1412.6980)を採用します．

```python
# 損失関数（DiceLoss）
loss = DiceLoss()

# 評価指標（IoU）
metrics = [
  IoU(threshold=0.5),    
]

# 最適化手法（Adam）
optimizer = torch.optim.Adam([dict(params=network.parameters())])
```

モデルの学習における振る舞いを定義します．上記で定義した損失関数，評価指標，最適化手法を設定します．

```python
# 学習用のエポック
train_epoch = TrainEpoch(
  network,
  loss = loss,
  metrics = metrics,
  optimizer = optimizer
)

# 評価用のエポック
eval_epoch = ValidEpoch(
  network,
  loss = loss,
  metrics = metrics
)
```

100エポックでモデルを学習します．各エポックで， **DiceLoss** と **IoU**が算出され，これを後述で可視化します．

```python
for i in range(100):
  train_logs = train_epoch.run(train_dataloader) # 学習
  print(train_logs)

  eval_logs = eval_epoch.run(eval_dataloader) # 評価
  print(eval_logs)
```

# 実験結果

**DiceLoss** と **IoU** の推移を次に示します．損失関数のDiceLossは，学習が進むと減少していくことが確認できます．また，評価指標のIoUは，学習が進むと増加していくことが確認できます．学習は順調に推移し，いずれも一定値に収束しており，過学習（オーバフィッティング）は確認できませんでした．

[![Image from Gyazo](https://i.gyazo.com/d5e86a2d678b994ace2a459e4d42a0bd.png)](https://gyazo.com/d5e86a2d678b994ace2a459e4d42a0bd)

[![Image from Gyazo](https://i.gyazo.com/75cec575e089b08dcf56c1a5a85055ac.png)](https://gyazo.com/75cec575e089b08dcf56c1a5a85055ac)

DiceとIoUの最大値を次の表にまとめます（Diceは1からDiceLossを引いた値）．国土地理院のサイトには，Dice（F値）が **0.8** 以上になったと記載されていました．今回の実験では，訓練データは **0.931** であり，0.8を超える結果となりました．一方，評価データは **0.600** であり，0.8には到達できませんでした．国土地理院の結果は，評価データに対するスコアだと思われ，大きな差があります．まだまだ工夫の余地がありそうです．

|  | Dice Score | IoU Score |
|:-:|:-:|:-:|
| 訓練データ | 0.931  | 0.872  |
| 評価データ | 0.600  |  0.799 |

次に学習したモデルを利用してマスク画像（道路位置の予測）を生成します（上から順に32.png，240.png，615.png）．左はオリジナル画像，中央は正解のマスク画像，右は予測したマスク画像です．いずれも道路のおおまかな形状や配置を予測できています．しかし，道路の接続が切断されてしまうなど，不完全な箇所も多く存在します．正確な情報を必要としない地物調査などの用途には，適用可能な精度と言えそうです．

[![Image from Gyazo](https://i.gyazo.com/92dfcd1a811d4b53548bd2990ef314d9.jpg)](https://gyazo.com/92dfcd1a811d4b53548bd2990ef314d9)
<!-- 32.png -->

[![Image from Gyazo](https://i.gyazo.com/8d094344a1e1fbbc5b2a123a9b38ce93.jpg)](https://gyazo.com/8d094344a1e1fbbc5b2a123a9b38ce93)
<!-- 240.png -->

[![Image from Gyazo](https://i.gyazo.com/0317c0ed5e06e5637495f2ecbba58811.jpg)](https://gyazo.com/0317c0ed5e06e5637495f2ecbba58811)
<!-- 615.png -->

# まとめ 

本記事では，国土地理院が公開した空中写真データセットを利用して，空中写真から道路を抽出するセマンティック・セグメンテーションを実装してみました．国土地理院のサイトに記載されている「Dice（F値）が0.8以上」には到達できませんでしたが，道路のおおまかな形状と配置の認識が可能であることを示すことができました．今後は，U-Netとは異なるモデル（[U-Net++](https://arxiv.org/abs/1807.10165)など）の適用を検討するとともに，セマンティック・セグメンテーションを活用した交通事故分析などの応用に取り組みたいと思っています．

# 出典

- 国土地理院（2022）：CNNによる道路抽出のための教師画像データ，国土地理院技術資料 H1-No.17．
- "Microsoft COCO: Common Objects in Context", Lin, Tsung-Yi et. al., arXiv, 2014
- "SegNet: A Deep Convolutional Encoder-Decoder Architecture for Image Segmentation", Badrinarayanan, Vijay et. al., arXiv, 2015
- "U-Net: Convolutional Networks for Biomedical Image Segmentation", Ronneberger, Olaf et. al., arXiv, 2015
- "Deep Residual Learning for Image Recognition", He, Kaiming et. al, arXiv, 2015
- "Adam: A Method for Stochastic Optimization", Kingma, Diederik P et. al., arXiv, 2014
- "UNet++: A Nested U-Net Architecture for Medical Image Segmentation", Zhou, Zongwei et. al., arXiv, 2018