1. 概要

以前，航空機や衛星画像の建物のセグメンテーションをkerasで実行しました．

衛星画像のSegmentation（セグメンテーション）により建物地図を作成する．

本記事では，流行っているPytorchを用いて，前回と同じ航空機の撮像画像から建物のセグメンテーションを行うモデルを構築します．例えば，予測モデルの出力は以下のようになりました．

画像の取得，前処理，PytorchによるDataの準備やモデルの構築，作成したモデルの検証など，Pytorchやセグメンテーションが初めての方を対象としたため，かなり細かく紹介しています．そのため，長文となりましたので，慣れている方はポイントだけ見てください．

ここで用いたコードはGithubにアップしましたので，ご興味のある方は試してみてください．Jupyter lab(notebook)で実行できます．
ご参考になれば幸いです．

環境

本記事の実装環境は以下となります．

OS:Ubuntu: 18.04LTS
GPU:GeoForce GTX1070

Python: 3.7
Pytorch: 1.1.0

2. セグメンテーションモデル

セグメンテーション（正確には，Semantic Segmentation)に関する記事は多数あります．例えば，以下の記事ではセグメンテーションのモデルで有名なU-Netが紹介されています．

U-NetでPascal VOC 2012の画像をSemantic Segmentationする (TensorFlow)

前回の記事と同様に，上記を参考にU-Netを構築することは学習としては非常に有益ですが，今回は複数のセグメンテーションを試すことができる以下のモジュールを使ってみました．

Segmentation Models
MIT license

このモジュールはPytorchをベースとした複数のSegmentationのモデルが準備されており，そのモデルには，Unet, Linknet, FPN, PSPNet, PAN, DeepLabV3があります．
また，それぞれのモデルの重み付けのパラメータも用意されており，転移学習やFine tuningができるため，比較的少ない学習データであっても高い精度のモデルが構築できることが期待されます．
さらに，このモジュールを用いた機械学習のコンペティションであるKaggleでの好成績の実績が紹介されており，その使い方についても学ぶことができます．

kaggle: Segmentation models

ここでは，セグメンテーションの代表モデルであるU-Netによる航空機の撮像画像より建物のセグメンテーションのモデルの構築例を紹介します．

3. 航空画像の取得

ここで用いる画像データは，前回と同様に以下のサイトより取得します．

Inria Aerial Image Labeling Dataset

ここには，それぞれ405 km²のアノテーション（学習用の建物地図）のある学習画像(Austin, Chicago,Kitsap Country, Western Tyrol, Viennaの５都市）と，アノテーションのないテスト画像になります．また，画像の分解能は30cmです．人工衛星で30cmの高い分解能をもつ人工衛星は，DigitalGlobeのWorldView-3があり，直下観測（センサが地球を直下にみたとき）で31cmの分解能となります．

航空機の撮像画像，およびアノテーション画像の例は以下となります．

航空機　撮像画像

アノテーション画像
Inriaが提供する画像データ例（Chicago)

WorldViewのサンプル画像は以下にありますので，ご興味がありましたらご覧ください．

Satellite Image Corporation: WorldView-3

Inriaのデータはタブの”Download”をクリックするとガイドがでてきますので，メールアドレス等を入力しダウンロードしてください．
全部で約20GBの画像データになります．

取得したデータや，モデルのスクリプト（コード）のディレクトリ（フォルダ）構成は以下となります．
ダウロードした画像は，imageフォルダにtrainおよびtestを置きます．
その画像から，後ほど説明します”前処理”にて学習用および検証用のデータを作成します．

以下は，今回のコードでの各ファイルやデータの位置構成になります．必要に応じて，それぞれの環境にてコードを編集ください．

4. セグメンテーションの実行

4.1 取得データの前処理

取得した撮像画像およびアノテーション画像を確認します．
最初にデータのファイルリストを作成し，それぞれの表示させます．

#モジュールのimport
from PIL import Image
import os
import sys
import numpy as np
from natsort import natsorted

import cv2
import matplotlib.pyplot as plt

import glob
#ダウンロードした画像のフィアルリストの作成
x_train_files = glob.glob('./images/*')
y_train_files = glob.glob('./gt/*')

#画像の読み込み
im_image = Image.open(x_train_files[0])

#画像をarrayに変換
im_list_image = np.asarray(im_image)
#貼り付け
plt.imshow(im_list_image)
plt.show()

#アノテーション画像の読み込み
im_gt = Image.open(y_train_files[0])
#画像をarrayに変換
im_list_gt = np.asarray(im_gt)
#貼り付け
plt.imshow(im_list_gt)
#表示
plt.show()

出力

アノテーション画像の画像構成を確認します．

print(im_list_gt)
print(len(im_list_gt))
print(im_list_gt.shape)

出力

[[255 255 255 ... 255 255 255]
 [255 255 255 ... 255 255 255]
 [255 255 255 ... 255 255 255]
 ...
 [  0   0   0 ...   0   0   0]
 [  0   0   0 ...   0   0   0]
 [  0   0   0 ...   0   0   0]]
5000
(5000, 5000)

これより，建物がある位置（pixel)は255，ないところは0となっていることがわかります．
また，その画像サイズは5000 x 5000pixelであることを確認しました．
次に，この画像を小さいサイズの画像に分割します．ここでは，各250ピクセルのサイズに分割しました．

#分割する画像サイズ，リサイズ画像のサイズ指定
height = 250
width = 250
img_size= len(im_list_gt)

#画像の分割処理関数
def ImgSplit(im):
    buff = []
    # 縦の分割枚数
    for h1 in range(int(img_size/height)):
        # 横の分割枚数
        for w1 in range(int(img_size/width)):
            w2 = w1 * height
            h2 = h1 * width
            #print(w2, h2, width + w2, height + h2)
            c = im.crop((w2, h2, width + w2, height + h2))
            buff.append(c)
    return buff

#撮像画像の分割
for i in range(len(x_train_files)):
# 画像の読み込み
    im=Image.open(x_train_files[i])
    file_name = os.path.splitext(os.path.basename(x_train_files[i]))[0]

    #画像の分割処理の実行
    hi=0
    for ig in ImgSplit(im):
        hi=hi+1
        # 保存先フォルダの指定
        ig.save("./train_images/"+ str(file_name) + '_' + str(i) +'_' + str(hi) +".png")

#アノテーション画像の分割
for i in range(len(y_train_files)):
# 画像の読み込み
    im=Image.open(y_train_files[i])
    file_name = os.path.splitext(os.path.basename(y_train_files[i]))[0]

    #画像の分割処理の実行
    hi=0
    for ig in ImgSplit(im):
        hi=hi+1
        # 保存先フォルダの指定
        ig.save("./train_gt/"+ str(file_name) + '_' + str(i) +'_' + str(hi) +".png")

画像の分割方法については，以前の記事をご参考にしてください．
分割したファイルは，元の画像ファイルのファイル名（都市名+番号）をヘッドにし，分割した順番をつけています．

次に，分割したデータの一部を検証用として利用します．
InriaにはTest用の画像が用意されていますが，それにはアノテーション画像がありません．そのため，構築したモデルを用いた予測結果との比較ができないため，ここでは学習用画像の一部を検証用として用います．そのための処理として，取得した画像のどれを検証用として用いるのか処理します．

DATA_DIR = './image/train/'

#分割した画像データのディレクトリを指定する．
x_train_dir = os.path.join(DATA_DIR, 'train_images')
y_train_dir = os.path.join(DATA_DIR, 'train_gt')

#分割した画像データのファイルリストを作成する．
x_train_files = glob.glob(x_train_dir +'/*')
y_train_files = glob.glob(y_train_dir +'/*')

#分割した画像ファイルのファイル数を確認する．
print(len(x_train_files))

次に，作成したファイルリストはその順番が適当なため，ソート処理により順序を固定します．ファイルリストのソートはnatsortを利用しました．
natsortをインストールされていない方は，'pip install natsort'を実行してインストールしてください．

#natsortをインストールしていない場合は以下を実行する．
#!pip install natsort

from natsort import natsorted

X_train_files = natsorted(x_train_files)
Y_train_files = natsorted(y_train_files)

ここで，ソートされたファイルリストを確認し，撮像画像のファイルとアノテーション画像のファイル順序が同じであることを確認します．

print(X_train_files[0])
print(Y_train_files[0])

出力

./image/train/train_images/chicago1_10_1.png
./image/train/train_gt/chicago1_10_1.png

次に，検証用画像として画像データのなかでAustinの撮像画像を用いるための処理を行います．ここでは，Austinとしましたが，何を指定されてもかまいません．

import shutil
import glob
import os

def move_glob(dst_path, pathname, recursive=True):
    for p in glob.glob(pathname, recursive=recursive):
        shutil.move(p, dst_path)

#Austinの撮像画像を検証用の画像データとして移動する．
move_glob('./image/train/val_images', './image/train/train_images/austin*.png')

#Austinのアノテーション画像を検証のアノテーション画像として移動する．
move_glob('./image/train/val_gt', './image/train/train_gt/austin*.png')

では，分割した画像を確認します．

#撮像画像の読み込み
im_image = Image.open(x_train_files[0])

#画像をarrayに変換
im_list_image = np.asarray(im_image)
#貼り付け
plt.imshow(im_list_image)
#表示
plt.show()

#アノテーション画像の読み込み
im_gt = Image.open(y_train_files[0])

#画像をarrayに変換
im_list_gt = np.asarray(im_gt)
#貼り付け
plt.imshow(im_list_gt)
#表示
plt.show()

アノテーション画像は，そのままでは建物の部分が255，建物でない部分が0となり，Pytorch用にデータをセットすると，Classの分類にエラーが生じます．今回は建物の有無の2値分類ですので，建物の部分を１，建物でない部分0とする処理をアノテーション画像に対して行います．

#学習用アノテーション画像の二値化処理
for i in range(len(y_train_files)):
    im = Image.open(y_train_files[i])
    im_list = np.asarray(im)
    X = im_list
    Y = np.where(X>1,1,X)
    pil_img = Image.fromarray(Y)
    pil_img.save(y_train_files[i])

#検証用アノテーション画像の二値化処理
for i in range(len(y_val_files)):
    im = Image.open(y_val_files[i])
    im_list = np.asarray(im)
    X = im_list
    Y = np.where(X>1,1,X)
    pil_img = Image.fromarray(Y)
    pil_img.save(y_val_files[i])

処理後のアノテーション画像を数値を確認します．

#画像の読み込み
im_gt = Image.open(y_train_files[0])
#画像をarrayに変換
im_list_gt = np.asarray(im_gt)

print(im_list_gt)

出力

[[1 1 1 ... 1 1 1]
 [1 1 1 ... 1 1 1]
 [1 1 1 ... 1 1 1]
 ...
 [1 1 1 ... 0 0 0]
 [1 1 1 ... 0 0 0]
 [1 1 1 ... 0 0 0]]

アノテーション画像が0,1で構成されていることを確認しました．
これで，データの前処理が終了です．

4.2 Datasetの作成

Pytorchでの画像処理のためのDatasetを作成します．
まずは，各モジュールをインストールします．

import os
os.environ['CUDA_VISIBLE_DEVICES'] = '0'

import numpy as np
import cv2
import matplotlib.pyplot as plt

Pytorchのモジュールをインストールし，GPUの動作を確認します．

import torch

if not torch.cuda.is_available():
  raise Exception("GPU not availalbe. CPU training will be too slow.")

print("device name", torch.cuda.get_device_name(0))

こちらを実行すると，使用しているGPUの情報が出力されます．
次に，PytorchのDatasetとDataloaderのモジュールをimportします．

from torch.utils.data import DataLoader
from torch.utils.data import Dataset as BaseDataset

Datasetは以下で作成します．
ここで，建物の有無の2値分類のため，そちらを設定します．

class Dataset(BaseDataset):
    """CamVid Dataset. Read images, apply augmentation and preprocessing transformations.

    Args:
        images_dir (str): path to images folder
        masks_dir (str): path to segmentation masks folder
        class_values (list): values of classes to extract from segmentation mask
        augmentation (albumentations.Compose): data transfromation pipeline 
            (e.g. flip, scale, etc.)
        preprocessing (albumentations.Compose): data preprocessing 
            (e.g. noralization, shape manipulation, etc.)

    """

    CLASSES = ['unlabelled', 'building'] #変更

    def __init__(
            self, 
            images_dir, 
            masks_dir, 
            classes=None, 
            augmentation=None, 
            preprocessing=None,
    ):
        self.ids = os.listdir(images_dir)
        self.images_fps = [os.path.join(images_dir, image_id) for image_id in self.ids]
        self.masks_fps = [os.path.join(masks_dir, image_id) for image_id in self.ids]

        # convert str names to class values on masks
        self.class_values = [self.CLASSES.index(cls.lower()) for cls in classes]

        self.augmentation = augmentation
        self.preprocessing = preprocessing

    def __getitem__(self, i):

        # read data
        image = cv2.imread(self.images_fps[i])
        image = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)
        mask = cv2.imread(self.masks_fps[i], 0)

        # extract certain classes from mask (e.g. cars)
        masks = [(mask == v) for v in self.class_values]
        mask = np.stack(masks, axis=-1).astype('float')

        # apply augmentations
        if self.augmentation:
            sample = self.augmentation(image=image, mask=mask)
            image, mask = sample['image'], sample['mask']

        # apply preprocessing
        if self.preprocessing:
            sample = self.preprocessing(image=image, mask=mask)
            image, mask = sample['image'], sample['mask']

        return image, mask

    def __len__(self):
        return len(self.ids)

Datasetに設定された画像データを以下を実行することで閲覧します．

def visualize(**images):
    """PLot images in one row."""
    n = len(images)
    plt.figure(figsize=(16, 5))
    for i, (name, image) in enumerate(images.items()):
        plt.subplot(1, n, i + 1)
        plt.xticks([])
        plt.yticks([])
        plt.title(' '.join(name.split('_')).title())
        plt.imshow(image)
    plt.show()

では，Datasetを設定し画像を閲覧します．

dataset = Dataset(x_train_dir, y_train_dir, classes=['building'])

image, mask = dataset[6] # get some sample
visualize(
    image=image, 
    building_mask=mask.squeeze(),
)

出力

航空機の撮像画像とアノテーション画像がセットで作成されていることが確認されました．

次に，学習用画像の水増しのためのAugmentationを実行します．
Augmentationはいくつかの方法はありますが，Segmentation Modelsで例として紹介されているalbumentationsを用いました．albumentationsについては以下の記事で詳しく紹介されていますので，こちらをご参考にしてください．

画像データ拡張ライブラリ ~ albumentations ~

まず，albumentationsをインストールします．pip install albumentationsでもよいそうですが,私の環境ではエラーが発生したため，以下の方法でインストールしました．エラーが発生する場合は試してみてください．

!pip install -U git+https://github.com/albu/albumentations --no-cache-dir

import albumentations as albu

def get_training_augmentation():
    train_transform = [

        albu.HorizontalFlip(p=0.5),

        albu.ShiftScaleRotate(scale_limit=0.5, rotate_limit=0, shift_limit=0.1, p=1, border_mode=0),

        albu.PadIfNeeded(min_height=320, min_width=320, always_apply=True, border_mode=0),
        albu.RandomCrop(height=320, width=320, always_apply=True),

        albu.IAAAdditiveGaussianNoise(p=0.2),
        albu.IAAPerspective(p=0.5),

        albu.OneOf(
            [
                albu.CLAHE(p=1),
                albu.RandomBrightness(p=1),
                albu.RandomGamma(p=1),
            ],
            p=0.9,
        ),

        albu.OneOf(
            [
                albu.IAASharpen(p=1),
                albu.Blur(blur_limit=3, p=1),
                albu.MotionBlur(blur_limit=3, p=1),
            ],
            p=0.9,
        ),

        albu.OneOf(
            [
                albu.RandomContrast(p=1),
                albu.HueSaturationValue(p=1),
            ],
            p=0.9,
        ),
    ]
    return albu.Compose(train_transform)


def get_validation_augmentation():
    """Add paddings to make image shape divisible by 32"""
    test_transform = [
        albu.PadIfNeeded(384, 480)
    ]
    return albu.Compose(test_transform)


def to_tensor(x, **kwargs):
    return x.transpose(2, 0, 1).astype('float32')


def get_preprocessing(preprocessing_fn):
    """Construct preprocessing transform

    Args:
        preprocessing_fn (callbale): data normalization function 
            (can be specific for each pretrained neural network)
    Return:
        transform: albumentations.Compose

    """

    _transform = [
        albu.Lambda(image=preprocessing_fn),
        albu.Lambda(image=to_tensor, mask=to_tensor),
    ]
    return albu.Compose(_transform)

多くの機能がありますが，まずは例で紹介されているAugmentationを実行しました．実行後の画像は以下で確認できます．

augmented_dataset = Dataset(
    x_train_dir, 
    y_train_dir, 
    augmentation=get_training_augmentation(), 
    classes=['building'],
)

# same image with different random transforms
for i in range(3):
    image, mask = augmented_dataset[8]
    visualize(image=image, mask=mask.squeeze(-1))

出力

変形された撮像画像に対して，アノテーション画像も同様に変形しているがわかります．
これで，Datasetの設定が終了です．

4.3 DataloaderとModelの作成

次に学習モデルと，モデルに入力するDataloaderを設定します．
今回は学習モデルにUnet,　学習済みモデルとしてResnet50を用いた転移学習を行います．

import torch
import numpy as np
import segmentation_models_pytorch as smp

ENCODER = 'resnet50'
ENCODER_WEIGHTS = 'imagenet'
CLASSES = ['building']
ACTIVATION = 'sigmoid' # could be None for logits or 'softmax2d' for multicalss segmentation
DEVICE = 'cuda'

# create segmentation model with pretrained encoder
model = smp.Unet(
    encoder_name=ENCODER, 
    encoder_weights=ENCODER_WEIGHTS, 
    classes=len(CLASSES), 
    activation=ACTIVATION,
)

preprocessing_fn = smp.encoders.get_preprocessing_fn(ENCODER, ENCODER_WEIGHTS)

Segmentation Modelsで利用できるモデルおよび学習済みモデルのリストは公式のGithubに記載されていますのでご参考ください．

次にDataloaderを設定します．

train_dataset = Dataset(
    x_train_dir, 
    y_train_dir, 
    augmentation=get_training_augmentation(), 
    preprocessing=get_preprocessing(preprocessing_fn),
    classes=CLASSES,
)

valid_dataset = Dataset(
    x_valid_dir, 
    y_valid_dir, 
    augmentation=get_validation_augmentation(), 
    preprocessing=get_preprocessing(preprocessing_fn),
    classes=CLASSES,
)

train_loader = DataLoader(train_dataset, batch_size=8, shuffle=True, num_workers=12)
valid_loader = DataLoader(valid_dataset, batch_size=1, shuffle=False, num_workers=4)

次は，損失係数です．Dice_lossを指標としてIoUも求めました．
Dice係数，IoUについては以下のサイトに詳しく説明されていますので，ご参考にしてください．

Deep Learning等の精度評価において、F値（Dice）とIoU（Jaccard）のどちらを選択するべきか？

loss = smp.utils.losses.DiceLoss()
metrics = [
    smp.utils.metrics.IoU(threshold=0.5),
]

活性化関数はAdamを用いました．

optimizer = torch.optim.Adam([ 
    dict(params=model.parameters(), lr=0.0001),
])

次に学習ループの設定です．今回は公式の紹介どおりとしています．

loss = smp.utils.losses.DiceLoss()
metrics = [
    smp.utils.metrics.IoU(threshold=0.5),
]

活性化関数は今回はAdamを用いました．

train_epoch = smp.utils.train.TrainEpoch(
    model, 
    loss=loss, 
    metrics=metrics, 
    optimizer=optimizer,
    device=DEVICE,
    verbose=True,
)

valid_epoch = smp.utils.train.ValidEpoch(
    model, 
    loss=loss, 
    metrics=metrics, 
    device=DEVICE,
    verbose=True,
)

最後に学習の実行です．以下にて学習をスタートします．
ここではエポックを40とし，25エポック時に学習レートを変更しています．
また，学習状況を可視化するために，各エポックでの学習および検証用データによるDice lossおよびIoUの記録を追加しています．

max_score = 0


#train accurascy, train loss, val_accuracy, val_loss をグラフ化できるように設定．
x_epoch_data = []
train_dice_loss = []
train_iou_score = []
valid_dice_loss = []
valid_iou_score = []


for i in range(0, 40):

    print('\nEpoch: {}'.format(i))
    train_logs = train_epoch.run(train_loader)
    valid_logs = valid_epoch.run(valid_loader)

    x_epoch_data.append(i)
    train_dice_loss.append(train_logs['dice_loss'])
    train_iou_score.append(train_logs['iou_score'])
    valid_dice_loss.append(valid_logs['dice_loss'])
    valid_iou_score.append(valid_logs['iou_score'])

    # do something (save model, change lr, etc.)
    if max_score < valid_logs['iou_score']:
        max_score = valid_logs['iou_score']
        torch.save(model, './best_model_Unet_resnet50.pth')
        print('Model saved!')

    if i == 25:
        optimizer.param_groups[0]['lr'] = 1e-5
        print('Decrease decoder learning rate to 1e-5!')

私の環境では，1エポックに45分ほどかかります．40エポックだと30時間ほどの計算時間となります．（2080tiがあれば...）

4.4 モデルの検証．

学習したモデルを検証する．
まず，学習時のDice_lossおよびIoUの推移を確認する．

fig = plt.figure(figsize=(14, 5))

ax1 = fig.add_subplot(1, 2, 1)
line1, = ax1.plot(x_epoch_data,train_dice_loss,label='train') 
line2, = ax1.plot(x_epoch_data,valid_dice_loss,label='validation')
ax1.set_title("dice loss")
ax1.set_xlabel('epoch')
ax1.set_ylabel('dice_loss')
ax1.legend(loc='upper right')

ax2 = fig.add_subplot(1, 2, 2)
line1, = ax2.plot(x_epoch_data,train_iou_score,label='train')
line2, = ax2.plot(x_epoch_data,valid_iou_score,label='validation') 
ax2.set_title("iou score")
ax2.set_xlabel('epoch')
ax2.set_ylabel('iou_score')
ax2.legend(loc='upper left')

plt.show()

出力．

学習時はエポック数に応じてDice損失は低下し，またIoUが向上しているのがわかる．一方，学習に用いていない検証画像の各指標は，よくはなっているが学習モデルほどではない．準備した学習画像の汎用性が低いことが疑われるため，Augmentationの追加や，学習用と検証用の画像を変えてみるなど，試してみるとわかるかもしれません．

学習したモデルを以下でロードし，モデルの確認を行います．

# load best saved checkpoint
best_model = torch.load('./best_model_Unet_resnet50.pth')

best_model

上記を実行すると，以下のようなモデルの詳細が閲覧できます．ここでは，一部のみ記載します．

Unet(
  (encoder): ResNetEncoder(
    (conv1): Conv2d(3, 64, kernel_size=(7, 7), stride=(2, 2), padding=(3, 3), bias=False)
    (bn1): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
    (relu): ReLU(inplace)
    (maxpool): MaxPool2d(kernel_size=3, stride=2, padding=1, dilation=1, ceil_mode=False)

では，構築したモデルでどこまで建物のセグメンテーションができているか，確認していきます．
公式サイトでは，検証用とはべつに試験用のデータを用いていますが，アノテーション画像がある試験用データがないため，ここでは検証用データを用います．

# create test dataset
test_dataset = Dataset(
    x_valid_dir, #x_test_dir
    y_valid_dir, #y_test_dir
    augmentation=get_validation_augmentation(), 
    preprocessing=get_preprocessing(preprocessing_fn),
    classes=CLASSES,
)

test_dataloader = DataLoader(test_dataset)

# evaluate model on test set
test_epoch = smp.utils.train.ValidEpoch(
    model=best_model,
    loss=loss,
    metrics=metrics,
    device=DEVICE,
)

logs = test_epoch.run(test_dataloader)

出力．

valid: 100%|██████████| 14400/14400 [05:30<00:00, 43.56it/s, dice_loss - 0.2963, iou_score - 0.6048]

最終結果を確認します．IoUが0.6程度なので，Inriaのサイトのテスト結果の基準指標と同程度ですね．これをベースに工夫することで，IoUの向上を図りたいと思います．IoUのベストスコアが0.8を超えているのは凄いです．

改めてDatasetを作成し，撮像画像，アノテーション画像および予測画像を比較します．

# test dataset without transformations for image visualization
test_dataset_vis2 = Dataset(
    x_valid_dir, y_valid_dir, #x_test_dir, y_test_dir,
    augmentation=get_validation_augmentation(), 
    classes=CLASSES,
)

for i in range(10):
    n = np.random.choice(len(test_dataset))

    image_vis = test_dataset_vis2[n][0].astype('uint8')
    image, gt_mask = test_dataset[n]

    gt_mask = gt_mask.squeeze()

    x_tensor = torch.from_numpy(image).to(DEVICE).unsqueeze(0)
    pr_mask = best_model.predict(x_tensor)
    pr_mask = (pr_mask.squeeze().cpu().numpy().round())

    visualize(
        image=image_vis, 
        ground_truth_mask=gt_mask, 
        predicted_mask=pr_mask
    )

出力．

結果みると，多くの建物は識別できているように見えますが，認識されていなかったり，もしくは誤検出している部分も見られます．これをみると，学習画像を工夫すればIoUを向上できそうです．

5. まとめ

Pytorchによる航空機の撮像画像から建物地図をつくるセグメンテーションの方法を紹介しました．

学習モデルをイチから作ることで理解が深まり，応用が広がりますが，まずは広く用いられているモデルを参考に簡易に構築できる方法として，Segmentation Modelsを利用しました．このモデルは設定できるパラメータが多いので，実験的にいろいろ試し，理解を深めるのによいツールと思います．

構造物のセグメンテーションやその地図であれば，Google Mapなどのサービスがすでにあります．私が航空機や人工衛星などの遠隔からの撮像画像を用いたかったのは，”変化”を捉えるためです．例えば，無料で提供されている欧州のSentinel-2の撮像は同一地点を5日毎に行われており，被雲率を考慮し月に１回利用したとしても，関心地域のトレンド（例えば街の発展具合など）を把握することができます．
また，災害前後の画像を比較すると，建物などの構造物の被害がフォーカスされることで，その被害規模を把握することに利用されることが期待されます．

次は，低解像度だけども広域を一度に撮像する人工衛星の撮像画像を対象に，セグメンテーションを実験していきます．

長文記事を最後までご覧いただきありがとうございました．私はこの分野の専門ではないため，間違って解釈しているところがあるかと思います．ご指摘いただければ幸いです．
また，コメント等ありましたらいただけると嬉しいです．

航空機や人工衛星の画像に関心を持つ方が多くなり，その応用の発展に寄与できれば幸いです．

6. 参考記事

衛星画像のSegmentation（セグメンテーション）により建物地図を作成する．
U-NetでPascal VOC 2012の画像をSemantic Segmentationする (TensorFlow)
Segmentation Models
MIT license
画像データ拡張ライブラリ ~ albumentations ~
Deep Learning等の精度評価において、F値（Dice）とIoU（Jaccard）のどちらを選択するべきか？

Pytorchによる航空画像の建物セグメンテーションの作成方法.