目次
今回は人工知能、ディープラーニングで細胞の種類を判別してみたいと思っています。判別するのは私が普段実験で使っている上皮細胞のEpH4という細胞と、遺伝子発現用のHEK293という細胞です。
EpH4
HEK293
見ると結構見た目が違います。この特徴を学習させます。
機械学習のフレームワークにはいろいろなものがありますが、今回はchainerを使います。インストールなどはLP-techの他の記事を参考にしてみるといいかもしれません。
この記事は畳み込みニューラルネットワークCNNを応用する内容です。
この記事は畳み込みニューラルネットワーク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の方が少し危なかったですが、どちらも正しく判定されています!!
終わりに
今回は機械学習で細胞の種類を判定する内容でした。これだけだと、やってみた、というだけの内容ですが、細胞の薬剤処理による形態変化の定量化などの様々な応用が考えられると思います。
機械学習はほとんどの特徴を学習できるのではないでしょうか。今後ますます活躍していく技術だと思っています。最後まで読んでいただき、ありがとうございました。
機械学習はほとんどの特徴を学習できるのではないでしょうか。今後ますます活躍していく技術だと思っています。最後まで読んでいただき、ありがとうございました。