衛星画像分析コンテストチュートリアル

create date : Jul. 23, 2018 at 20:37:23

これは「産業技術総合研究所　衛星画像分析コンテスト」( https://signate.jp/competitions/61 )のチュートリアルです。分析する際の参考になれば幸いです。

データをダウンロードして深層学習モデリングによって投稿ファイルを出力するところまでできるようになることを期待しています。このチュートリアルの残りの構成は以下のようになります。

  1. データのダウンロードとファイルの解凍
  2. データの読み込みと可視化
  3. 画像データの前処理
  4. モデリング
  5. 予測結果の出力
  6. まとめ

想定している分析環境は以下です。

  • OS: Windows 10 Pro
  • Anaconda==5.0.1(Python==3.6.3)
  • Library
    • chainer==3.3.0
    • cupy==2.3.0 (chainerをGPU上で動かすために必要)
    • matplotlib==2.1.0
    • numpy==1.13.3
    • pandas==0.20.3
    • scikit-image==0.13.0
    • tqdm==4.19.5
    • その他の依存ライブラリ
    • GPU: Quadro M1200 (optional)

1. データのダウンロードとファイルの解凍

  1. https://signate.jp/howto に従い、会員登録を行う。
  2. https://signate.jp/competitions/61/data へアクセスし、データ一式（train_1.zip, train_2.zip, train_3.zip, test_1.zip, test_2.zip, test_3.zip, test_4.zip, train_master.tsv）をダウンロードする。
  3. ダウンロードしたファイルを作業ディレクトリにおいて、zipファイルを解凍する。すると学習用画像データフォルダtrainと評価用画像データフォルダtestが作成される。

2. データの読み込みと可視化

まずは作業ディレクトリに移動し、画像データ(train_1.tif)を読み込んでみます。読み込むにはskimageを用います。

In [1]:
import os
from skimage import io
In [2]:
image_path = os.path.join('train', 'train_1.tif')
image = io.imread(image_path)
In [3]:
print(image.shape)

(32, 32, 7)

データ構造は(縦, 横, チャンネル)という形になっています。https://signate.jp/competitions/61 にあるように、今回のデータの場合は縦×横=32×32で、チャンネル数は7となっていることが確認されました。

次に実際に画像データを可視化してみます。それぞれのチャンネルについて可視化します。例えば１番目のチャンネルのデータを取り出したいときはimage[:,:,1]のように記述します。可視化にはmatplotlibを用います。

In [4]:
from matplotlib import pyplot as plt
%matplotlib inline
In [5]:
fig, (ax0, ax1, ax2, ax3, ax4, ax5, ax6) = plt.subplots(nrows=1,ncols=7,figsize=(10, 3))
ax0.imshow(image[:,:,0])
ax0.set_title('1')
ax0.axis('off')
ax0.set_adjustable('box-forced')

ax1.imshow(image[:,:,1])
ax1.set_title('B')
ax1.axis('off')
ax1.set_adjustable('box-forced')

ax2.imshow(image[:,:,2])
ax2.set_title('G')
ax2.axis('off')
ax2.set_adjustable('box-forced')

ax3.imshow(image[:,:,3])
ax3.set_title('R')
ax3.axis('off')
ax3.set_adjustable('box-forced')

ax4.imshow(image[:,:,4])
ax4.set_title('5')
ax4.axis('off')
ax4.set_adjustable('box-forced')

ax5.imshow(image[:,:,5])
ax5.set_title('6')
ax5.axis('off')
ax5.set_adjustable('box-forced')

ax6.imshow(image[:,:,6])
ax6.set_title('7')
ax6.axis('off')
ax6.set_adjustable('box-forced')

fig.tight_layout()
In [6]:
import pandas as pd
data = pd.read_csv('train_master.tsv', sep='\t')
data.head()
Out[6]:
file_name flag
0 train_0.tif 0
1 train_1.tif 0
2 train_2.tif 0
3 train_3.tif 0
4 train_4.tif 0

train_1.tifのflagは0なので、読み込んだtrain_1.tifはゴルフ場を全く含まない画像データでした。

3. 画像データの前処理

モデリングを行う際、多くの場合、学習時または予測時に画像データに対して前処理を行います。今回は画像データとして格納されている値の範囲を0と1の間に収まるような前処理を行います。そのためにskimageのexposureモジュールのrescale_intensityを用います。このメソッドは読み込んだ画像データを渡すと値の範囲が0と1の間となるデータを返します。詳しくは http://scikit-image.org/docs/stable/api/skimage.exposure.html#skimage.exposure.rescale_intensity をご覧ください。

In [7]:
from skimage import exposure
In [8]:
image_rescaled = exposure.rescale_intensity(image)
In [9]:
# 前処理を行う前
print('最大値:', image.max())
print('最小値:', image.min())

最大値: 23325.0
最小値: 6652.0
In [10]:
# 前処理を行った後
print('最大値:', image_rescaled.max())
print('最小値:', image_rescaled.min())

最大値: 1.0
最小値: 0.0

4. モデリング

畳み込みニューラルネットワークモデルを構築する

画像データに対してゴルフ場が含まれているか否かを判定する畳み込みニューラルネットワークモデルを構築します。今回はオープンソースの深層学習フレームワーク の一つであるchainer( https://chainer.org/ )を用います。今回実装するのは前半が畳み込み層2つ、後半が全結合層2つからなる比較的浅いネットワークです。実装についてさらに詳しく知りたい人はhttps://docs.chainer.org/en/stable/tutorial/ をご覧ください。

In [11]:
import chainer.functions as F
import chainer.links as L
from chainer import Chain
In [12]:
class Cnn(Chain):
    def __init__(self, out_size):
        super(Cnn, self).__init__(
                conv1 = L.Convolution2D(None, 32, 3, stride=2),
                conv2 = L.Convolution2D(None, 32, 3, stride=2),
                fc3 = L.Linear(None, 512),
                fc4 = L.Linear(None, out_size))

        self.train = True
    
    def __call__(self, x):
        with chainer.using_config('enable_backprop', self.train):
            h = F.max_pooling_2d(F.relu(self.conv1(x)), 2, stride=2)
            h = F.max_pooling_2d(F.relu(self.conv2(h)), 2, stride=2)
            h = F.relu(self.fc3(h))
            y = self.fc4(h)
        
            return y

今回は、「ゴルフ場が含まれている」か「ゴルフ場が含まれていない」かの2クラス分類モデルを構築します。よって、out_size=2とします。

In [13]:
model = Cnn(out_size=2)

モデルを学習させる

処理の進捗確認するためのライブラリtqdmをインポートします。

In [14]:
from tqdm import tqdm_notebook as tqdm
tqdm.monitor_interval = 0

1枚の画像データを渡して画像処理後のデータを出力する関数を実装します。

In [15]:
def preprocess(image, mode = 'train'):
    """
    image: shape = (h, w, channel)を想定。
    mode: 'train', 'val', 'test'を想定。
    """
    if mode == 'train':
        # その他いろいろな前処理メソッドを実装してみてください
        if image.max()!=image.min():
            image = exposure.rescale_intensity(image)

    elif mode == 'val':
        # その他いろいろな前処理メソッドを実装してみてください
        if image.max()!=image.min():
            image = exposure.rescale_intensity(image)

    elif mode == 'test':
        # その他いろいろな前処理メソッドを実装してみてください
        if image.max()!=image.min():
            image = exposure.rescale_intensity(image)
    else:
        # その他いろいろな前処理メソッドを実装してみてください
        if image.max()!=image.min():
            image = exposure.rescale_intensity(image)
        
    return image

ミニバッチ数分のデータを読み込み、前処理されたデータを返す関数を実装します。ラベルが存在する場合（学習時'train'、検証時'val'）には画像データとともにラベルも返します。

In [16]:
def generate_minibatch(data_path, minibatch_meta, mode = 'train'):
    images = []
    if mode == 'train' or mode=='val':
        labels = []
    for data in minibatch_meta.iterrows():
        im_path = os.path.join(data_path, data[1]['file_name'])
        image = io.imread(im_path)

        # preprocess image
        image = preprocess(image, mode = mode)
        image = image.transpose((2,0,1))
        
        if mode == 'train' or mode=='val':
            labels.append(data[1]['flag'])

        images.append(image)

    images = np.array(images)
    if mode == 'train' or mode=='val':
        labels = np.array(labels)
        
        return images, labels
    else:
        return images

データを学習用と検証用に分ける関数を実装します。ratioで学習用の比率を指定します。デフォルトでは95%を学習用として分割します。乱数生成などにnumpyを使うのでインポートします。

In [17]:
import numpy as np
In [18]:
def split_data(data, ratio=0.95):
    train_index = np.random.choice(data.index, int(len(data)*ratio), replace=False)
    val_index = list(set(data.index).difference(set(train_index)))
    train = data.iloc[train_index].copy()
    val = data.iloc[val_index].copy()
    
    return train, val

今回の評価関数となっているIOU(intersection over union)を実装します。正解データy_trueと予測されたデータy_predを渡してIOUを返します。

In [19]:
def IOU(y_true, y_pred):
    y_true = np.array(y_true)
    y_pred = np.array(y_pred)
    p_true_index = np.where(y_true==1)[0]
    p_pred_index = np.where(y_pred==1)[0]
    union = set(p_true_index).union(set(p_pred_index))
    intersection = set(p_true_index).intersection(set(p_pred_index))
    if len(union)==0:
        return 0
    else:
        return len(intersection)/len(union)

ミニバッチサイズを指定して、その枚数分画像を読み込み、前処理を行い、モデルの学習パラメータを更新していく関数を実装します。デフォルトで一回につき64枚処理するようにします(batchsizeを参照)。損失関数としてはsoftmax cross entropy( https://en.wikipedia.org/wiki/Cross_entropy )を用います。今回は正例と負例の比率が偏っているため、損失関数を計算する際に偏りを解消するためにclass_weightを設定します。class_weightは正例と負例の比率の逆比とします(F.softmax_cross_entropyに渡すclass_weightを参照)。過剰適合を防ぐためにL2正則化を行います。係数は0.001とします(chainer.optimizer.WeightDecayを参照)。最適化手法としてはNesterovの方法( https://arxiv.org/abs/1212.0901 )を適用します。学習率は0.01とします(init_lrを参照)。ゴルフ場が含まれるかどうかの判断に用いる閾値はゴルフ場が含まれるクラスの確率を小さい順にソートし、後ろから数えて正例の数番目の確率値としています(th_trとth_valを参照)。GPUを使う場合はgpu=Trueとします(GPU環境がない場合はgpu=Falseとします)。

In [20]:
import chainer
import chainer.cuda as cuda
from chainer import optimizers, serializers
In [21]:
def train_model(model, train, val, data_path, model_name, batchsize = 64, init_lr = 0.01, num_epochs = 30, gpu = True):
    optimizer = optimizers.NesterovAG(lr = init_lr) # 学習率はinit_lr
    optimizer.setup(model)
    optimizer.add_hook(chainer.optimizer.WeightDecay(0.001)) # L2正則化
    
    best_score = 0
    best_th = 0.5
    th_tr = 0.5
    th_val = 0.5
    tr_num_positives = train['flag'].sum()
    val_num_positives = val['flag'].sum()

    if gpu:
        model.to_gpu()
    epbar = tqdm(range(num_epochs))
    for epoch in epbar:
        epbar.set_description('epoch {}'.format(epoch+1))
        # training
        batch_index = np.random.permutation(np.arange(len(train))//batchsize)
        train['batch_index'] = batch_index
        num_samples = 0
        count = 0
        train_loss = 0
        train_iou = 0
        tr_proba = np.array([])
        tr_true = np.array([])
        pbar = tqdm(train.groupby('batch_index'), leave=False)
        pbar.set_description('training')
        for t in pbar:
            model.cleargrads()
            minibatch, labels = generate_minibatch(data_path, t[1], mode = 'train')
            weight = np.array([sum(labels)/len(labels), 1-sum(labels)/len(labels)]).astype(np.float32)
            minibatch = minibatch.astype(np.float32)
            labels = labels.astype(np.int32)
            if gpu:
                minibatch = cuda.to_gpu(minibatch)
                labels = cuda.to_gpu(labels)
                weight = cuda.to_gpu(weight)
            
            y = model(minibatch)
            loss = F.softmax_cross_entropy(y, labels, class_weight = weight)
            loss.backward()
            optimizer.update()
            y_proba = cuda.to_cpu(F.softmax(y).data)[:,-1]
            y_pred = (y_proba > th_tr).astype(np.int)
            tr_proba = np.concatenate((tr_proba, y_proba))
            labels = cuda.to_cpu(labels)
            tr_true = np.concatenate((tr_true, labels))
            train_loss += float(loss.data)*len(minibatch)
            train_iou += IOU(labels, y_pred)
            num_samples += len(minibatch)
            count += 1
            pbar.set_postfix(IOU=round(train_iou/count,4), loss=round(train_loss/num_samples, 4), th=round(th_tr, 2))
        
        pbar.close()
        train_iou = IOU(tr_true, (tr_proba>th_tr).astype(np.int))
        
        # validation
        val_proba = np.array([])
        val_true = np.array([])
        batch_index = np.random.permutation(np.arange(len(val))//batchsize)
        val['batch_index']=batch_index
        count = 0
        val_iou = 0
        model.train = False
        pbar = tqdm(val.groupby('batch_index'), leave=False)
        pbar.set_description('validating')
        for v in pbar:
            minibatch, labels = generate_minibatch(data_path, v[1], mode = 'val')
            minibatch = minibatch.astype(np.float32)
            labels = labels.astype(np.int32)
            if gpu:
                minibatch = cuda.to_gpu(minibatch)
            y = model(minibatch)
            y_proba = cuda.to_cpu(F.softmax(y).data)[:,-1]
            y_pred = (y_proba > th_val).astype(np.int)
            val_proba = np.concatenate((val_proba, y_proba))
            val_true = np.concatenate((val_true, labels))
            count+=1
            val_iou += IOU(labels, y_pred)
            pbar.set_postfix(IOU=round(val_iou/count, 4), th=round(th_val, 2))
        pbar.close()
        val_iou = IOU(val_true, (val_proba>th_val).astype(np.int))
        
        
        tr_proba.sort()
        th_tr = tr_proba[-tr_num_positives]
        val_proba.sort()
        th_val = val_proba[-val_num_positives]

        if best_score <= val_iou:
            best_score = val_iou
            best_th = th_val
            serializers.save_npz(model_name.split('.')[0]+'_best.npz', model)
        serializers.save_npz(model_name, model)
        epbar.write('[epoch {}] train_iou:{}  val_iou:{}  best_iou:{}  best_th:{}'.format(epoch+1, round(train_iou,4), round(val_iou,4),
                                                                                          round(best_score,4), round(best_th,4)))
        
        model.train = True

    print('best IOU:', best_score, 'best threshold:', best_th)
    return model

実際にモデルを学習させてみます。まず、作業ディレクトリに移動し、ダウンロードしたデータを読み込みます。次に学習用データを学習用と検証用に改めて分割します。比率は学習用データを95%とします。そしてモデルを学習させます。エポック数(学習データ全体を学習させる繰り返し回数)は20とします(num_epochsを参照)。

In [22]:
# 学習用画像が格納されているディレクトリを指定する
data_path = 'train'

# 学習用データを学習用と検証用に改めて分割する
train, val = split_data(data, ratio = 0.95)

print('-'*20, 'train', '-'*20)
print('number of samples:', len(train))
print('number of positives:', train['flag'].sum())
print('nubmer of negatives:', (1- train['flag']).sum())
print('-'*47)

print('-'*20, 'val', '-'*20)
print('number of samples:', len(val))
print('number of positives:', val['flag'].sum())
print('nubmer of negatives:', (1- val['flag']).sum())
print('-'*45)

-------------------- train --------------------
number of samples: 281372
number of positives: 11793
nubmer of negatives: 269579
-----------------------------------------------
-------------------- val --------------------
number of samples: 14810
number of positives: 649
nubmer of negatives: 14161
---------------------------------------------
In [ ]:
# モデルを学習させる
model_name = 'Cnn.npz'
model = train_model(model, train, val, data_path, model_name, num_epochs=20)

想定している環境では約3時間16分かかります。

5. 予測結果の出力

学習済みモデルと評価用画像データを渡して予測結果を返す関数を実装します。

In [23]:
def predict(model, data_path, mode = 'test', batchsize = 128, gpu = True, threshold = 0.5):
    data = pd.DataFrame()
    data['file_name'] = os.listdir(data_path)
    predictions = np.array([])
    model.train = False
    if gpu:
        model.to_gpu()
    data['batch']=np.arange(len(data))//batchsize
    for t in tqdm(data.groupby('batch')):
        minibatch = generate_minibatch(data_path, t[1], mode = mode)
        minibatch = minibatch.astype(np.float32)
        if gpu:
            minibatch = cuda.to_gpu(minibatch)
        y = model(minibatch)
        y_proba = cuda.to_cpu(F.softmax(y).data)[:,-1]
        y_pred = (y_proba>threshold).astype(np.int)
        predictions = np.concatenate((predictions, y_pred))
    data['prediction'] = predictions
    del data['batch']
    
    return data

学習済みモデルを用いて評価用データに対して予測結果を出力します。まず評価用データが格納されているディレクトリを指定します。ゴルフ場が含まれているかの判断のための閾値thresholdは検証の結果、大きめが良さそうなので0.9にしてみます。提出の形式はtsvでヘッダーなしなので、そのように設定します。GPUを使う場合はgpu=Trueとします(GPU環境がない場合はgpu=Falseとします)。

In [ ]:
# 評価用データが格納されているディレクトリを指定し、閾値を0.9に設定する
data_path = 'test'
threshold = 0.9

# 学習済みのモデルをロードする
model = Cnn(out_size=2)
model_name = 'Cnn_best.npz'
serializers.load_npz(model_name, model)

# 予測結果を出力する
submit = predict(model, data_path, threshold=threshold)

想定している環境では約18分かかります。

In [ ]:
submit[['file_name', 'prediction']].to_csv('submit.tsv', sep='\t', header=None, index=None)

投稿すると、スコアは0.1～0.2ほどになります。

6. まとめ

以上「産業技術総合研究所　衛星画像分析コンテスト」で深層学習モデリングを行うための実装について一通り説明しました。まだまだ工夫の余地はあります。例えば(1)リアルタイムで学習データを水増ししてバリエーションを増やす、(2)予測時に評価用データに対して複数種類の前処理を行い、それぞれに対する結果を平均化する(test time augmentation)、(3)より性能の良い深層学習モデル(ResNetなど)を利用する、(4)確率値に対する閾値を最適化する、(5)異なるモデルをアンサンブルするなどがあげられます。まだコンペティション終了まで時間はあるのでこのチュートリアルを参考にしつつ工夫をしてみてください。皆様のご応募をお待ちしております。

create date : Jul. 23, 2018 at 20:37:23