[実践編]Deep Learningによる衛星画像の建物Segmentation.

先日，衛星画像から建物のSegmentationの方法，およびその対象画像の前処理，後処理に関する記事を書きました．
・衛星画像のSegmentation（セグメンテーション）により建物地図を作成する．
・Deep Learning(深層学習)のための画像の標準化，リサイズ，分割，結合処理．
今回は，これらのコードを総合して，対象画像から建物のセグメンテーションまで通して処理できるようにしました．
また，テスト用としてコードおよび学習モデルも公開します．よろしければ，試してみてください．

##1.　環境設定．
主な環境は以下の通りで，これにて対象とする衛星画像（テスト画像）の建物Segmentationを行っています．GPUがないマシンでも，その機能を使わなければ，時間はかかりますが処理は可能と思います．

環境
・ubuntu 16.04LTS
・Cuda 8.0
・GPU GTX1070

・Python 3.5
・keras 2.1.5
・numpy 1.13.3
・natsort5.5.0
・tensorflow-gpu 1.5.0
・tensorflow-gpu-base 1.5.0
・tensorflow-tensorboad 1.5.1

##2. 衛星画像の建物Segmentation．
　前回は，インプットデータ，教師データを用いた学習モデルの作成から行いました．今回は，これまで学習したモデルデータを用いて，対象の衛星画像（テストデータ）から，建物セグメンテーションの結果のみを作成します．
　以下がメインのコードになります．

main.py

import os
from matplotlib import pyplot
from keras.optimizers import Adam

from unet import UNet
from dice_coefficient import dice_coef_loss, dice_coef
from images_loader import load_images, save_images
from option_parser import get_option

#バックエンドをインポート
from keras.backend import tensorflow_backend as backend

from PIL import Image
import sys
import numpy as np
from natsort import natsorted

#分割する画像サイズ，リサイズ画像のサイズ指定
height = 256
width = 256
H_img_size=8960
V_img_size=6144

INPUT_IMAGE_SIZE = 256
BATCH_SIZE = 10
EPOCHS = 1000

DIR_MODEL = os.path.join('..', 'Model')
DIR_OUTPUTS = os.path.join('..', 'OutputsA')
DIR_TESTS = os.path.join('..', 'TestDataA')

File_MODEL = 'detect_satellite_model13.hdf5'


def predict(input_dir):
    (file_names, inputs) = load_images(input_dir, INPUT_IMAGE_SIZE)

    network = UNet(INPUT_IMAGE_SIZE)
    model = network.model()
    model.load_weights(os.path.join(DIR_MODEL, File_MODEL))
    preds = model.predict(inputs, BATCH_SIZE)

    save_images(DIR_OUTPUTS, preds, file_names)

#画像の規格化処理関数
def normalize(arr):
    arr = arr.astype('float')
    # Do not touch the alpha channel
    for i in range(3):
        minval = arr[...,i].min()
        maxval = arr[...,i].max()
        if minval != maxval:
            arr[...,i] -= minval
            arr[...,i] *= (255.0/(maxval-minval))
    return arr

#画像の分割処理関数
def ImgSplit(im):
    # 読み込んだ画像を２５６*２５６のサイズに分割する

    buff = []
    # 縦の分割枚数
    for h1 in range(int(V_img_size/height)):
        # 横の分割枚数
        for w1 in range(int(H_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

#画像の結合処理関数
def ImageMerge():
    image_list = os.listdir(DIR_OUTPUTS)
    #print(image_list)
    image_list_sort1=sorted(image_list)
    #print(image_list_sort1)
    image_list_sort=natsorted(image_list)
    #print(image_list_sort)

    nh=int(H_img_size/width)
    nv=int(V_img_size/height)

    total_width = width * nh
    max_height = height * nv

    new_im = Image.new("RGB", (int(total_width), int(max_height)))

    h = 0
    w = 0

    for i in range(int(nh*nv)):
        new_im.paste(Image.open(DIR_OUTPUTS + "/"+ image_list_sort[i]), (w*width,h*height))
        w = w +1
        if w % int(nh) == 0:
           w = w-int(nh)
           h = h +1

    new_im.save(DIR_OUTPUTS + '/test_total.jpg')

if __name__ == '__main__':
    #画像の読み込み
    im=Image.open(DIR_TESTS+"/test.TIF")

    #画像のリサイズ・標準処理
    img_resize=im.resize((H_img_size,V_img_size))
    arr = np.array(img_resize)
    new_img = Image.fromarray(normalize(arr).astype('uint8'))

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

    args = get_option(EPOCHS)
    EPOCHS = args.epoch

    if not(os.path.exists(DIR_MODEL)):
        os.mkdir(DIR_MODEL)
    if not(os.path.exists(DIR_OUTPUTS)):
        os.mkdir(DIR_OUTPUTS)

    #画像ファイルの解析実行　(DIR_INPUTS)
    predict(DIR_TESTS)

    #画像結合処理の実行
    ImageMerge()

    #解析結果は，OutPutsAフォルダのtotal_test.jpgとして保存される．

#処理終了時に下記をコール
backend.clear_session()

対象とする衛星画像を，TestDataAフォルダにいれmain.pyコードを実行します．
ここでは，その衛星画像(test.TIF)のサイズが9137☓6391サイズであったため,できるだけそのサイズのまま使いたいので，幅と高さを256ピクセルの整数倍である8960☓6144にリサイズして用いました．これで，840に分割されてSegmentation処理が行われます．
結果は，OutPutsAフォルダにtotal_test.jpgとして出力されます．参考に．途中の840の画像ファイルも同じく保存しています．
　サンプルの衛星画像，学習モデルも含めた一連のコードファイルをこちらに置きました．よろしければ，ダウロードし実行してみてください．
衛星画像の建物Segmentationファイル一式

##3. 衛星画像の建物Segmentation例．
　今回は，学習用に用いた航空写真画像ではなく，高分解能の衛星画像をサンプル画像として用いました．
　高分解能の衛星画像の価格はかなり高いため，フランスの光学観測衛星SPOTのサンプル画像を使っています．
衛星画像のサンプル画像は，以下のAirBusのホームページよりダウンロードできます．
[AIRBUS Defense and Space sample imagery] (https://www.intelligence-airbusds.com/en/8262-sample-imagery)
0.5mの分解能のPleiadesの画像もありましたが，まずは1.5mの分解能のSPOTの画像で解析してみました．
　ページの右の”Availabel Products”からOptical Imageryを選択し，次に衛星名のSPOTを選択します．
　複数のサンプル衛星画像がありますが，ここでは”SPOT 7 Ortho Display Pansharpened (GEOTIFF 8bits)”を選び，イタリアのベニスの観測画像をサンプル画像として解析してみました．
　サンプル画像，およびその建物Segmentationの結果は以下となります．

サンプル衛星画像（ベニス）

建物セグメンテーション解析結果

建物はある程度は識別できてそうかな，というレベルです．細かく見てみます．

対象画像

建物セグメンテーション解析結果

こうみると，それほど建物が識別できているように見えません．特に，影がある部分は検出できていないのがわかります．
学習モデルの画像は，雲や建物の影がない理想的なものでした．それに対して，今回　対象とした衛星画像は一般的な衛星画像であり，その場合はこのように検出率（精度）が低い結果となっています．
　これより，学習モデルに用いる教師画像は，対象とする衛星画像により作成した，できるだけ多くの条件を含む画像が必要である，ということがわかりました．これより，解析精度（今回は，建物の検出率）は，Segmentation手法やコードそのものよりも，どういった学習画像を用いてモデルをつくるか，に大きく影響することがわかりました．なので，よく”教師データが重要”と言われるのですね．
　衛星画像の教師データはなかなか入手できないため，自分で作成するしかないかと思います．サンプルデータはある程度ありますので，今後はアノテーションソフトを用いて，地道に教師データ（建物のマスク画像）を作ってみようと思います．

##4.　まとめ．
　今回は，航空写真画像の教師データにより作成した学習モデルを用いて，衛星画像の規格化，リサイズ，分割の前処理，解析画像の結合処理まで，ひとつのコードで処理してみました．よろしければ，みなさんも一度試してみてください．ご参考になれば幸いです．