MNISTのニューラルネットワークを使ってDICOM画像の部位を分類してみる
Keyword：DICOM，Chainer，MNIST
3.01でMNISTのネットワークを用いて手書き文字の分類を行いました。
3.02ではMNISTの手書き文字のデータ構造を解析してDICOM画像の訓練データの作り方を調べました。
そして3.03でDICOM画像を訓練データの画像リストファイルととテストデータの画像リストファイルに分けるプログラムを作成しました。
この章では，これらのプログラムを使って訓練データファイルからCT画像の部位（Head，Lung，Chest，Abdomenの4種類）の分類問題を学習，テストするプログラムを作成します。
実験用のDICOM画像はここにあります。私のCT画像です。実験にのみ使用し，他の目的で使用しないでください。
プログラムのポイント
- 前提として./TrainTestListディレクトリに訓練データファイルTrainData.0.3.csvとテストデータファイルTestData.0.3.csvが存在するとします。
- この章では学習プログラムdicomNNtrain.pyとテストプログラムdicomNNtest.pyを作成します。これらのプログラムは前章のmnistNN.py(Pyゼミ3.01参照)のネットワークモデルクラスを使用しますが，さらに分類するためのクラスClassifyImageクラスを作成しました。
- 学習プログラムdicomNNtrain.pyは訓練データファイルTrainData.0.3.csvの中に記述された画像ファイルを用いて学習を行い，学習が終了すると，学習モデルをMyModel.hdf5にファイル出力します。
- テストプログラムdicomNNtest.pyは学習モデルMyModel.hdf5ファイルを読み込み，TestData.0.3.csvファイルに指定されたテストDICOM画像について分類テストを行います。
学習プログラム（dicomNNtrain.py）の概要
import
DICOM画像を読み込むためpydicomが追加されています。constractTrainData関数
この関数はPyゼミ3.02のプログラムと同じです。訓練用のDICOM画像のリストファイル（画像ファイル名，ラベル番号とラベル名のリスト）を引数に受け取り，訓練データを作成して返します。
readDicom2png関数
この関数はPyゼミ1.02のプログラムと同じです。DICOM画像を適正な濃度にウィンドニングして画素値を返します。
getLabelName関数
Pyゼミ3.03で作成したラベル名のCSVファイルlabelName.csvを読み込み，ラベル数（分類すべき数）と配列（ラベル番号とラベル名）を返します。MyModelクラス
ニューラルネットワークを定義するクラスです。
Pyゼミ3.01の手書き文字認識のMyModelクラスとほとんど同じですが，init関数がラベル数nlblを引数に持ち，出力層l3に設定しています。
手書き文字認識は10種類の数字ですが，今回のCT画像の部位分類では4種類になります。
l3 = L.Linear( 100, nlbl )
call関数はPyゼミ3.01と同じです。
ClassifyImageクラス
今回学習のために新しく作成したクラスです。２つの関数を持ちます。
init関数はニューラルネットワークモデルの初期化と設定を行っています。
ClassifierはChainerが提供する関数で，引数に与えたネットワークモデルに対して損失や精度を計算します。
ここでは損失関数にソフトマックスークロスエントロピイF.softmax_cross_entropyを指定しています。
doTrain関数は訓練を実行するときに呼ばれる関数です。
trainer.run関数により訓練が行われ，終了後に学習モデルをファイルmodelNameに保存します。
chainer.serializers.save_hdf5( self.modelName, self.model )
プログラム実行の流れを見てみる
main関数を見てみましょう。
変数trainfnmとlabelfnmにそれぞれ訓練画像ファイルのリストファイル名とラベル名のファイルが代入されている。
画像のリサイズresizeの値を読み込みます。
CT画像は512×512ですが，MNISTのネットワークは28x28（入力ノード数784）なのでこの値を設定して画像をリサイズします。
学習後に保存するモデルの名前をmodelnmに設定します。
学習単位のデータ数であるbatch_sizeと全データの学習の回数n_epochをそれぞれ10と25に設定します。
次にgetLabelName関数を呼んで，ラベル数nLblとラベル名配列lblNameを得ます。
次にconstractTrainData関数を呼んで訓練データをtrainに代入します。
次に保存するモデル名とラベル数をClassifyImageクラスに与えてインスタンスciを作成します。
最後にClassifyImageクラスのdoTrain関数を呼び出して，訓練を行います。
プログラムを入力して実行してみよう
プログラムは約110行ありますが，これまでのPyゼミ3.01から3.03で扱ってきた関数です。
新しいコードは少ないです。
dicomNNtrain.py
import os, sys
import numpy as np
import cv2
import chainer
from chainer import Function, report, training, utils, Variable
from chainer import datasets, iterators, optimizers, serializers
from chainer import Link, Chain, ChainList
import chainer.functions as F
import chainer.links as L
from chainer.datasets import tuple_dataset
from chainer.training import extensions
import pydicom
def constractTrainData(fname , resz):
image,label = [], []
for i, line in enumerate(open(fname, 'r')):
data = line.split(",") # 画像ファイル名と正解データに分割
print("Train",i, data[0], "Label#",data[1])
img = readDicom2png(data[0], "tmp.png")
if resz > 0: #画像のリサイズ
img = cv2.resize(img,(resz, resz))
img = img.astype(np.float32)
img = (img - 128)/128 #画素値0-255を-1.0から+1.0に変換
image.append([img])
t = np.array(int(data[1]), dtype=np.int32)
label.append(t)
train = tuple_dataset.TupleDataset(image, label)
return train
def readDicom2png(dcmfnm, tmpfnm):
ds = pydicom.read_file(dcmfnm)
wc = ds.WindowCenter
ww = ds.WindowWidth
img = ds.pixel_array
max = wc + ww / 2
min = wc - ww / 2
img = 255 * (img - min)/(max - min)
img[img > 255] = 255
img[img < 0 ] = 0
img = img.astype(np.uint8)
cv2.imwrite(tmpfnm, img)
return img
def getLabelName(fname):
lblnm = []
for line in open(fname, 'r'):
data = line.split(",")
lblnm.append(data[1].strip())
print("LabelName:",lblnm)
return len(lblnm), lblnm
class MyModel(Chain): #ネットワークモデルを定義する
def __init__(self, nlbl):
super(MyModel, self).__init__(
l1 = L.Linear(784,100),
l2 = L.Linear(100,100),
l3 = L.Linear(100,nlbl),
)
def __call__(self, x): #順伝播計算
h1 = F.relu(self.l1(x)) #活性化関数にReLuを定義
h2 = F.relu(self.l2(h1))
return self.l3(h2)
class ClassifyImage():
def __init__(self, mdlnm, nlbl):
model = MyModel(nlbl)
self.model = L.Classifier(model,lossfun=F.softmax_cross_entropy)
self.opt = optimizers.Adam()
self.opt.setup(self.model)
self.modelName = mdlnm
def doTrain(self, nLbl, train, batch_size, n_epoch, modelName):
iterator = iterators.SerialIterator(train, batch_size, shuffle=True)
updater = training.StandardUpdater(iterator, self.opt)
trainer = training.Trainer(updater, (n_epoch, 'epoch'), out='result')
trainer.extend(extensions.LogReport())
trainer.extend(extensions.PrintReport(['epoch', 'main/loss', 'main/accuracy']))
trainer.extend(extensions.ProgressBar())
trainer.run()
print("学習モデル",self.modelName,"を保存します．.. . .")
chainer.serializers.save_hdf5( self.modelName, self.model )
print("\t... . .学習モデルを保存しました")
if __name__=='__main__':
trainfnm = "./TrainTestList/TrainData.0.3.csv"
labelfnm = "./TrainTestList/labelName.csv"
resize = 28 #画像リサイズ
modelnm = "MyModel.0.3.hdf5" #モデル保存名
batch_size = 10 # バッチサイズ
n_epoch = 15 # エポック数
nLbl, lblName = getLabelName(labelfnm)
train = constractTrainData(trainfnm, resize)
print("\nTrain data information\n")
print("画像データ数:",len(train))
print("画像/ラベル：",len(train[0]),"\n画像Chanel数:",len(train[0][0]))
print("高さ:",len(train[0][0][0]),"\n幅 :",len(train[0][0][0][0]))
ci = ClassifyImage(modelnm, nLbl)
ci.doTrain(nLbl, train, batch_size, n_epoch, modelnm)
実行結果
最初にラベル名labelName.csvを読み込みラベル名をターミナルに表示します。
その後訓練画像を読み込み，画像ファイル名とラベル番号をターミナルに表示します。
訓練データ情報としてデータ数142や画像のサイズ28ｘ28などが表示され，訓練が始まります。
学習を重ねるうちに損失（main/loss）が低下し，精度（main/accuracy）が向上していることがわかります。
最終的に学習モデルMyModel.0.3.hdf5が保存されます。
$ python dicomNNtrain.py
LabelName: ['Lung', 'Abdomen', 'Chest', 'Head']
Train 0 ../dcmdir2/Lung/W5305890 Label# 0
Train 1 ../dcmdir2/Lung/K5303593 Label# 0
Train 2 ../dcmdir2/Lung/R5305640 Label# 0
Train 3 ../dcmdir2/Lung/T5304000 Label# 0
Train 4 ../dcmdir2/Lung/A5304343 Label# 0
Train 5 ../dcmdir2/Lung/D5304468 Label# 0
・・・・・
Train 137 ../dcmdir2/Head/D5221046 Label# 3
Train 138 ../dcmdir2/Head/F5221140 Label# 3
Train 139 ../dcmdir2/Head/G5221187 Label# 3
Train 140 ../dcmdir2/Head/Q5221671 Label# 3
Train 141 ../dcmdir2/Head/E5221093 Label# 3
Train data information
画像データ数: 142
画像/ラベル： 2
画像Chanel数: 1
高さ: 28
幅 : 28
epoch main/loss main/accuracy
1 0.807113 0.713333
2 0.222487 0.95
3 0.106435 0.992857
4 0.07357 0.985714
5 0.0498589 0.992857
・・・・・
21 0.0171901 0.993333
22 0.00569069 1
23 0.0104065 0.992857
24 0.0069965 0.992857
25 0.0046676 1
学習モデル MyModel.0.3.hdf5 を保存します．.. . .
... . .学習モデルを保存しました
テストプログラム（dicomNNtest.py）の作成
訓練プログラムで保存された学習モデルのファイルを読み込み，さらにテスト画像を読み込み部位の分類のテストを行うプログラムを作成します。
テストプログラムの概要
import
訓練プログラムと同様にpydicomをインポートします。constructTestData関数
テスト画像ファイルのリストファイル（画像ファイル名，ラベル番号とラベル名のリスト）から画像データ,ラベル番号と画像ファイル名の配列testを返します。imgは2次元配列の画像，data[1]はラベル番号を整数に変換しています。
data[0]は画像ファイル名です。
この3つを配列の要素として，配列testに追加しています。
test.append( [ img, int( data[1] ), data[0] ] )
readDicom2png関数
訓練プログラムdicomNNtrain.pyの関数と同じです。getLabelName関数
訓練プログラムdicomNNtrain.pyの関数と同じです。MyModelクラス
訓練プログラムdicomNNtrain.pyの関数と同じです。ClassifyImageクラス
init関数でモデルを定義し，学習モデルmyModel.hdf5を読み込みます。predict関数はテスト画像testimgを引数に持ち，np.arrayに変換したあと，ChainerのVariable型に変換します。
predictorメソッドを使って，推論を行いSoftmax関数を用いて確率に変換します。
img = np.array( [ testimg ] )
x = Variable( img )
y = self.model.predictor( x )
y = F.softmax( y )
確率のベクトルyで最も高い値をargmaxメソッド使って求めて推論結果（ラベル番号）返します。
answer = y.data
answer = np.argmax(answer, axis=1)
プログラム実行の流れを見てみる
main関数を見てみましょう。
変数testfnmとlabelfnmにそれぞれテスト画像ファイルのリストファイル名とラベルファイル名を代入します。
resizeはMNISTの手書き画像にリサイズするため28に設定します。
modelnmは訓練プログラムで保存されたモデル名を設定します。
次にgetLabelName関数を呼んで，ラベル数をnLblにラベル名を配列lblNameに格納します。
配列lblNameのインデックスはラベル番号に一致する。つまりラベル番号0のラベル名はlblName[0]になります。
各ラベルの正答数を計数するため配列を作成します。
例えば，ラベル番号が1のデータが正答するとlblCnt[1]をインクリメントします。
例えば，ラベル番号が1のデータが正答するとlblCnt[1]をインクリメントします。
lblCnt = np.zeros( nLbl ) #各ラベルの正答数を計数
lblAll = np.zeros( nLbl ) #各ラベルの総数を計数
次にconstructTestData関数を呼び，すべてのテスト画像，正解ラベルと画像ファイル名を配列dataに格納します。
ClassifyImageクラスを呼びインスタンスciを作成します。
全体の正答数を計数する変数okを初期化します。
次にfor文にて配列dataからひとつづつデータを取り出し，テスト画像，正解ラベルと画像ファイル名をtestimg, ans, fnmに代入します。
テスト画像testimgをモデルciのpredict関数に渡し，推論を行い，推論結果predictionと確率分布softmaxを得ます。
もし推論結果predictionと正答ラベルansが一致していれば，その正答ラベルの配列をインクリメントします。
最後に，全体の平均正答率と個々のラベルの正答率を表示します。
プログラムを入力して実行してみよう
プログラムは100行を超えていますが，訓練プログラムと共通のところが多くあります。共通部分と異なる部分を意識しながら入力してみましょう。
import os, sys
import numpy as np
import cv2
import chainer
import chainer.links as L
import chainer.functions as F
from chainer import optimizers, Chain, Variable, initializers
from chainer.training import extensions
import pydicom
def constructTestData(fname , resz):
test = []
for i, line in enumerate(open(fname, 'r')):
data = line.split(",")
print("Test",i, data[0], "Label#",data[1])
img = readDicom2png(data[0], "tmp.png")
if resz > 0:
img = cv2.resize(img,(resz, resz))
img = img.astype(np.float32)
img = (img - 128)/128
test.append([[img], int(data[1]), data[0]])
return test
def readDicom2png(dcmfnm, tmpfnm):
ds = pydicom.read_file(dcmfnm)
wc = ds.WindowCenter
ww = ds.WindowWidth
img = ds.pixel_array
max = wc + ww / 2
min = wc - ww / 2
img = 255 * (img - min)/(max - min)
img[img > 255] = 255
img[img < 0 ] = 0
img = img.astype(np.uint8)
cv2.imwrite(tmpfnm, img)
return img
def getLabelName(fname):
lblnm = []
for line in open(fname, 'r'):
data = line.split(",")
lblnm.append(data[1].strip())
print("LabelName:",lblnm)
return len(lblnm), lblnm
class MyModel(Chain): #ネットワークモデルを定義する
def __init__(self, nlbl):
super(MyModel, self).__init__(
l1 = L.Linear(784,100),
l2 = L.Linear(100,100),
l3 = L.Linear(100,nlbl),
)
def __call__(self, x): #順伝播計算
h1 = F.relu(self.l1(x)) #活性化関数にReLuを定義
h2 = F.relu(self.l2(h1))
return self.l3(h2)
class ClassifyImage:
def __init__(self, modelname, nlbl):
model = MyModel(nlbl)
self.model = L.Classifier(model)
print("modelファイル ",modelname," を読み込みます．．．")
chainer.serializers.load_hdf5( modelname, self.model )
print("．．．ファイルを読み込みました")
def predict(self, testimg):
img = np.array([testimg])
#print(">img>>",type(img),len(img),len(img[0]),len(img[0][0]),len(img[0][0][0]))
x = Variable(img)
y = self.model.predictor(x)
y = F.softmax(y)
answer = y.data
answer = np.argmax(answer, axis=1)
return answer[0], y[0].data
if __name__=='__main__':
testfnm = "./TrainTestList/TestData.0.3.csv"
labelfnm = "./TrainTestList/labelName.csv"
resize = 28 #画像リサイズ
modelnm = "MyModel.0.3.hdf5" #モデル保存名
nLbl, lblName = getLabelName(labelfnm)
lblCnt = np.zeros(nLbl) #各ラベルの正答数を計数
lblAll = np.zeros(nLbl) #各ラベルの総数を計数
data = constructTestData(testfnm, resize)
print("data:",len(data),len(data[0]))
print("img :",len(data[0][0]), len(data[0][0][0]),len(data[0][0][0][0]))
ci = ClassifyImage(modelnm, nLbl) #モデルの読み込み
ok = 0
np.set_printoptions(formatter={'float': '{: 0.5f}'.format})
for i, line in enumerate(data):
testimg, ans, fnm = line
prediction, softmax = ci.predict(testimg) #推論
print("\n#",i,"Softmax:",softmax)
result = "File:"+fnm+"\tAns:"+lblName[ans]+"→ Predict:"+lblName[prediction]
lblAll[ans] += 1 #ラベルがansの数の計数
if prediction == ans: #推論結果と答えが一致したら
print(result + "\tOK")
ok += 1
lblCnt[prediction] += 1
else:
print(result + "\tNG")
print("\n平均正答率：",ok/len(data))
for i in range(nLbl):
print(str(i), lblName[i],"\t",lblCnt[i]/lblAll[i])
実行結果
最初にラベル名labelName.csvを読み込みラベル名をターミナルに表示します。
次にテスト画像と正解ラベル番号が読み込み，ターミナルに表示します。
その後，学習モデルを読み込み，各画像について推論を行い，正答ラベルと推論結果を表示します。
両者が一致していればOKが，そうでなければNGが表示されます。
最終的に平均正答率と各ラベルの正答率を表示します。
$ python dicomNNtest.py LabelName: ['Lung', 'Abdomen', 'Chest', 'Head'] Test 0 ../dcmdir2/Lung/E5306265 Label# 0 Test 1 ../dcmdir2/Lung/M5305390 Label# 0 Test 2 ../dcmdir2/Lung/C5304421 Label# 0 Test 3 ../dcmdir2/Lung/O5303765 Label# 0 Test 4 ../dcmdir2/Lung/P5305546 Label# 0 Test 5 ../dcmdir2/Lung/Q5305593 Label# 0 ・・・・・ Test 68 ../dcmdir2/Head/B5220953 Label# 3 Test 69 ../dcmdir2/Head/P5221625 Label# 3 Test 70 ../dcmdir2/Head/L5221421 Label# 3 Test 71 ../dcmdir2/Head/W5222046 Label# 3 Test 72 ../dcmdir2/Head/X5222093 Label# 3 data: 73 3 28 modelファイル myModel.0.3.hdf5 を読み込みます．．． ．．．ファイルを読み込みました # 0 Softmax: [ 0.99994 0.00006 0.00000 0.00000] File:../dcmdir2/Lung/E5306265 Ans:Lung→ Predict:Lung OK # 1 Softmax: [ 0.99971 0.00001 0.00024 0.00004] File:../dcmdir2/Lung/M5305390 Ans:Lung→ Predict:Lung OK # 2 Softmax: [ 0.99986 0.00000 0.00011 0.00002] File:../dcmdir2/Lung/C5304421 Ans:Lung→ Predict:Lung OK # 3 Softmax: [ 0.99996 0.00000 0.00003 0.00001] File:../dcmdir2/Lung/O5303765 Ans:Lung→ Predict:Lung OK ・・・・・ # 70 Softmax: [ 0.00002 0.00011 0.00002 0.99985] File:../dcmdir2/Head/L5221421 Ans:Head→ Predict:Head OK # 71 Softmax: [ 0.00023 0.00061 0.00009 0.99908] File:../dcmdir2/Head/W5222046 Ans:Head→ Predict:Head OK # 72 Softmax: [ 0.00084 0.00582 0.00023 0.99311] File:../dcmdir2/Head/X5222093 Ans:Head→ Predict:Head OK 平均正答率： 0.9863013698630136 0 Lung 1.0 1 Abdomen 1.0 2 Chest 0.933333333333 3 Head 1.0
このプログラムは訓練プログラムと共通する部分が多いのでプログラム概要は割愛します。
おわりに
mnistのニューラルネットワークをDICOM画像に発展して適用しました。
かなりシンプルなネットワークで，しかも本来512x512のCT画像を28x28にリサイズして学習していますが，正答率は良い結果が得られています。
（ただし，同一人物のCT画像なので上手く行って当然ですね）
次に，Deep Neural Networkの代表であるAlexnetについてPyゼミ3.05で実験してみたいと思います。今回のプログラムが完成すると，次はモデルを変えるだけです。
0 件のコメント:
コメントを投稿