細胞種を機械学習で判別する！ - IMACEL Academy -人工知能・画像解析の技術応用に向けて-| エルピクセル株式会社

機械学習のフレームワークにはいろいろなものがありますが、今回はchainerを使います。インストールなどはLP-techの他の記事を参考にしてみるといいかもしれません。

この記事は畳み込みニューラルネットワークCNNを応用する内容です。

今回は、まず、この2つの細胞の画像のうち右側の高密度の画像をもとに5＊5＝25枚の分割画像を作ります。次に各25枚すなわち、50枚の画像を3:1で学習用と検証用に分けます。画像数がやや少ないですがデモということで！（途中で画像の数を水増しする処理も紹介しています！！！）

まずはフォルダ構成から。

フォルダ構成

chainer.pyがpythonコードの書いてあるファイルです。

original_dataには元のデータを入れます。その中の0のフォルダにEpHの画像を、1のフォルダにHEK293の画像を入れています。今回は1枚ずつしかありませんが、何枚入れてもOKです。

画像の分割と明るさの統一

このコードが分割を行う関数になります。

import cv2
import os.path
import glob

#1枚の画像をN*Nの画像に分割、縮小もする、色補正
def img_split(path,N,label):
    files = glob.glob(path)
    files.sort()
    for i, filename in enumerate(files):
        #print filename, os.path.isfile(filename)
        print i,
        img = cv2.imread(filename, cv2.IMREAD_GRAYSCALE)
        img = cv2.resize(img,(64*N,64*N))
        img = img * 64.0 / cv2.mean(img)[0]
        height, width = img.shape
        
        for x in range(0,N):
            for y in range(0,N):
                clp = img[height/N*y:height/N*(y+1), width/N*x:width/N*(x+1)]
                output = "../train_data/"+label+"/"+str(i*N*N+x*N+y)+'.png'
                cv2.imwrite(output, clp)

分割と明るさの統一.py

img_split(r"../original_data/0/*.png",5,"0")
img_split(r"../original_data/1/*.png",5,"1")

分割.py

画像の水増し

分割しただけでは学習に十分な数の画像とは言えないため、回転によって数を増やします。

import cv2
import numpy as np

def rotateImage(image, angle):
    center = tuple(np.array(image.shape)[:2]/2) # 2d !
    rot_mat = cv2.getRotationMatrix2D(center,angle,1.0)
    result = cv2.warpAffine(image, rot_mat, tuple(np.array(image.shape)[:2]) ,flags=cv2.INTER_LINEAR)
    return result

def rotateImageAndSave(path, label):
    files = glob.glob(path)
    files.sort()
    for i, filename in enumerate(files):
        print i,
        img = cv2.imread(filename, cv2.IMREAD_GRAYSCALE)
        for deg in range(90,360,90):
            img = rotateImage(img, deg)
            str_deg = str(deg)
            output = "../train_data/"+label+"/"+str(i)+"rotate"+str_deg+'.png'
            cv2.imwrite(output, img)
            flipX = cv2.flip(img, 0)
            output = "../train_data/"+label+"/"+str(i)+"rotate"+str_deg+"flipX"+'.png'
            cv2.imwrite(output, flipX)
            flipY = cv2.flip(img, 1)
            output = "../train_data/"+label+"/"+str(i)+"rotate"+str_deg+"flipY"+'.png'
            cv2.imwrite(output, flipY)
            flipXY = cv2.flip(img, -1)
            output = "../train_data/"+label+"/"+str(i)+"rotate"+str_deg+"flipXY"+'.png'
            cv2.imwrite(output, flipXY)

水増し関数.py

rotateImageAndSave(r"../train_data/1/*.png","1")
rotateImageAndSave(r"../train_data/0/*.png","0")

水増し.py

画像パスとラベルの対応関係を記述する

inputfolder = r"../train_data"

import sys
import os.path
import shutil
import glob
import random

# 書き込むテキストファイル
train_list = open("train.txt", "w")
test_list = open("test.txt", "w")

classno = 0
count = 0
for label in range(0,2):
    label=str(label)
    imagefolder = inputfolder+"/"+label+"/"

    images = [r.split('/')[-1] for r in glob.glob(imagefolder+"*.png")]
    startcount = count
    length = len(images)
    for image in images:
        # 各クラスのデータの75%を学習に25%をテストに使う
        if random.random() < 0.75:
            train_list.write(imagefolder+image+" %d\n" %classno)
        else:
            test_list.write(imagefolder+image+" %d\n" %classno)
        count+=1
    classno+=1

train_list.close()
test_list.close()
print "Complete!"

ラベル対応.py

学習モデル

from chainer import cuda, Function, FunctionSet, gradient_check, Variable, optimizers, serializers
from chainer.training import extensions
import chainer.functions as F

class CNN():
    #cuda.get_device(0).use()    # GPU用コード
    def __init__(self):
        self.model = FunctionSet(
            conv1 = F.Convolution2D(1,32,3,pad=1),
            conv2 = F.Convolution2D(32,32,3,pad=1),
            conv3 = F.Convolution2D(32,16,3,pad=1),
            conv4 = F.Convolution2D(16,16,3,pad=1),
            l1 = F.Linear(256,50),#1つ目の値は計算方法があるが適当でもChainerが実行時に正しい値を教えてくれる
            l2 = F.Linear(50,3),
        ).to_cpu()
            

    def forward(self, x_data, y_data, train=True):
        x, t = Variable(x_data), Variable(y_data)
        c = F.max_pooling_2d(F.relu(self.model.conv1(x)), 2)
        c = F.max_pooling_2d(F.relu(self.model.conv2(c)), 2)
        c = F.max_pooling_2d(F.relu(self.model.conv3(c)), 2)
        c = F.max_pooling_2d(F.relu(self.model.conv4(c)), 2)
        h = F.dropout(F.relu(self.model.l1(c)), train=train)
        y = self.model.l2(h)

        if train:   # 学習時
            return F.softmax_cross_entropy(y, t)
        else:
            return y       # 評価時

モデル.py

いよいよ学習！

学習の前に元データの画像サイズを確認しておきます。

今回は64と出力されました。1片64の正方形画像です。

import cv2
img = cv2.imread(r"../train_data/0/0.png", 0)
print img.shape
imagesize = img.shape[0]

画像サイズ.py

学習のためのコードです。今回は100回学習を繰り返します。

import numpy as np
import chainer
from chainer import Function, FunctionSet, gradient_check, Variable, optimizers, serializers
from chainer.training import extensions
import argparse
import random
import chainer.functions as F
import cv2
import time

#cuda.get_device(0).use()    # GPU用コード

# 学習モデル設定
model = CNN()
optimizer = optimizers.Adam()
optimizer.setup(model.model)

# 学習データリストファイルから一行ずつ読み込む
train_list = []
for line in open(r"train.txt"):
    pair = line.strip().split()
    train_list.append((pair[0], np.float32(pair[1])))

# 画像データとラベルデータを取得する
x_train = []    # 画像データ格納
y_train = []    # ラベルデータ格納
for filepath, label in train_list:
    img = cv2.imread(filepath, 0)   # グレースケールで読み込む
    x_train.append(img)
    y_train.append(label)

# 学習で使用するsoftmax_cross_entropyは
# 学習データはfloat32,ラベルはint32にする必要がある。
x_train = np.array(x_train).astype(np.float32)
y_train = np.array(y_train).astype(np.int32)
# 画像を（学習枚数、チャンネル数、高さ、幅）の4次元に変換する
x_train = x_train.reshape(len(x_train), 1, imagesize, imagesize) / 255

N = len(y_train)

batchsize = 60
datasize = len(x_train)
n_epoch =  100   # 繰り返し学習回数
# 学習開始
train_start = time.time()

print "start"
for epoch in range(1, n_epoch+1):
    sum_loss = 0
    perm = np.random.permutation(N) # データセットの順番をシャッフル
    train_start_one = time.time()
    for i in range(0,datasize, batchsize):
        x_batch = x_train[perm[i:i+batchsize]] # バッチサイズ分のデータを取り出す
        y_batch = y_train[perm[i:i+batchsize]]
        # 初期化
        optimizer.zero_grads()
        # 誤差伝播
        loss = model.forward(x_batch, y_batch, train=True)
        # 誤差逆伝播
        loss.backward()
        # ネットワークパラメータ更新
        optimizer.update()

        sum_loss += loss.data
    print "epoch:{} loss={} time:{}".format(epoch, sum_loss / (datasize/batchsize), time.time()-train_start_one)
    if epoch%10==0:
        serializers.save_npz("model_{}".format(n_epoch), model.model)
        serializers.save_npz("state_{}".format(n_epoch), optimizer)
        
print "train time:{}".format(time.time()-train_start)
# 学習したモデルを出力する
serializers.save_npz("model_{}".format(n_epoch), model.model)
serializers.save_npz("state_{}".format(n_epoch), optimizer)

学習.py

学習完了！

学習過程

私のMacbook Air ではtrain time:1088.80467081とでました。約15分くらいですね。続いて精度の検証を行ってみたいと思います。

精度の検証

学習に使っていない25％の画像を使って検証を行います。

import matplotlib.pyplot as plt

def draw_digit3(data, n, ans, recog):
    plt.subplot(20, 20, n)
    Z = data.reshape(imagesize,imagesize)
    Z = Z[::-1,:]             # flip vertical
    plt.xlim(0,imagesize-1)
    plt.ylim(0,imagesize-1)
    plt.pcolor(Z)
    if ans==recog:
        plt.title("ans=%d, recog=%d"%(ans,recog), size=8, color="c")
    else:
        plt.title("ans=%d, recog=%d"%(ans,recog), size=8, color="r")
    plt.gray()
    plt.tick_params(labelbottom="off")
    plt.tick_params(labelleft="off")
    
plt.figure(figsize=(30,30))

import numpy as np
import chainer
from chainer import cuda, Function, FunctionSet, gradient_check, Variable, optimizers, serializers
from chainer.training import extensions
import argparse
import random
import chainer.functions as F
import cv2
import time

#cuda.get_device(0).use()    # GPU
print "test start"

# 学習モデル読み込み
model = CNN()
serializers.load_npz("model_100", model.model)

# テストデータリストファイルから一行ずつ読み込む（学習時と同じ）
test_list = []
for line in open(r"../train_data/index/test.txt"):
    pair = line.strip().split()
    test_list.append((pair[0], np.float32(pair[1])))

# 画像データとラベルデータを取得する（学習時と同じ）
x_test = []    # 画像データ格納
y_test = []    # ラベルデータ格納
for filepath, label in test_list:
    img = cv2.imread(filepath, 0)   # グレースケールで読み込む
    x_test.append(img)
    #print filepath
    y_test.append(label)

# 学習で使用するsoftmax_cross_entropyは
# 学習データはfloat32,ラベルはint32にする必要がある。
x_test = np.array(x_test).astype(np.float32)
y_test = np.array(y_test).astype(np.int32)
# 画像を（学習枚数、チャンネル数、高さ、幅）の4次元に変換する
x_test = x_test.reshape(len(x_test), 1, imagesize, imagesize) / 255

N = len(y_test)

batchsize = 1
datasize = len(x_test)

# 判定開始
test_start = time.time()
false_count = 0
perm = np.random.permutation(N) # データセットの順番をシャッフル
for i in range(0,datasize, batchsize):
    x_batch = x_test[[i]] # バッチサイズ分のデータを取り出す
    y_batch = y_test[[i]]

    result = model.forward(x_batch, y_batch, train=False) # 学習ではないのでtrainをFalse
    if (i%2==0):
        draw_digit3(x_batch[0], (i)/2+1, y_batch[0], np.argmax(result.data))    
    if np.argmax(result.data) != y_batch[0]:
        false_count += 1
        print "{} No.{} NG! correct:{}, result:{}".format(test_list[i][0], i, y_batch[0], np.argmax(result.data))
    else:
        print "{} No.{} OK! correct:{}, result:{}".format(test_list[i][0], i, y_batch[0], np.argmax(result.data))

print "data num:{} false num:{} accuracy={}".format(datasize, false_count, 1 - (float(false_count) / datasize))
print "test time:{}".format(time.time()-test_start)

検証.py

100%の精度で判別できています。

精度検証

まあ、元が同じ画像だからね。。

というあなたのために、低密度の方の画像をEpH4かHEK293か判定してみましょう。

低密度EpH4

低密度HEK293

import numpy as np
import chainer
from chainer import cuda, Function, FunctionSet, gradient_check, Variable, optimizers, serializers
from chainer.training import extensions
import argparse
import random
import chainer.functions as F
import cv2
import time
import os.path
import glob


# 学習モデル読み込み
model = CNN()
serializers.load_npz("model_100", model.model)

def judge(path,N):
    img_files = glob.glob(path)
    for i, filename in enumerate(img_files):
        img = cv2.imread(filename, cv2.IMREAD_GRAYSCALE)
        #print img.shape
        img = cv2.resize(img,(64*N,64*N))
        img = img * 64.0 / cv2.mean(img)[0]
        height, width = img.shape
        x_test = []
        for x in range(0,N):
            for y in range(0,N):
                clp = img[height/N*y:height/N*(y+1), width/N*x:width/N*(x+1)]
                x_test.append(clp)
        #print x_test[0].shape
        
        x_test = np.array(x_test).astype(np.float32)
        x_test = x_test.reshape(len(x_test), 1, 64, 64) / 255
        batchsize = 1
        datasize = len(x_test)
        # 判定開始
        test_start = time.time()
        true_count = 0;
        false_count = 0;
        for i in range(0,datasize):
            x_batch = x_test[[i]]
            y_batch = y_test[[i]]

            result = model.forward(x_batch, y_batch, train=False) # 学習ではないのでtrainをFalse
            if (np.argmax(result.data)==0):
                false_count +=1    
            else:
                true_count += 1
        print filename, int(100.0*float(true_count)/(true_count+false_count)),"%\t","(" ,true_count, "vs", false_count, ")\t",
        if true_count>false_count:
            print "HEK293"
        else:
            print "EpH4"

画像分類！.py

さて結果は？？

judge("../not_included/*.png",5)

ジャッジ！.py

判定結果

EpH4の方が少し危なかったですが、どちらも正しく判定されています！！

終わりに

今回は機械学習で細胞の種類を判定する内容でした。これだけだと、やってみた、というだけの内容ですが、細胞の薬剤処理による形態変化の定量化などの様々な応用が考えられると思います。

機械学習はほとんどの特徴を学習できるのではないでしょうか。今後ますます活躍していく技術だと思っています。最後まで読んでいただき、ありがとうございました。

参考文献

畳み込みニューラルネットワークについての詳しい解説

http://postd.cc/how-do-convolutional-neural-networks-work/