Pyゼミ3.03 DICOM画像を訓練用とテスト用に分けてみる

DICOM画像を訓練用とテスト用に分けてみる

Keyword:python, DICOM， クロスバリデーション


AIの実験では，ある訓練データの集合に対して学習を行って，その学習モデルを使って未知のテストデータに対して推論を行い，その正答率から学習モデルの性能を評価します。

この実験のためにたくさん集めたDICOM画像を訓練用とテスト用に用手的に振り分けるのは大変な作業になります。
そこで，ここではラベル付されたDICOM画像を訓練用とテスト用に分けるプログラムを作成します。
つまり，画像ファイル名とラベル番号のリストになった訓練用CSVファイルとテスト用CSVファイルを作成するプログラムを作成します。

準備：ラベル付されたディレクトリの中にDICOM画像を格納する

このプログラムは，
「振り分けるDICOM画像は，親ディレクトリの下にラベル付されたサブディレクトリがあり，その中に該当する複数の画像が保存されている。」
こととします。
サンプルのDICOM画像はここにあります。私の全身のCT画像です。個人の実験にのみ使用し，他では使わないでください。

このサンプルDICOM画像では，部位を分類するタスクを考えます。
親ディレクトリdcmdir2の下に，Abdomen（腹部），Chest（胸部），Head（頭部）とLung（肺野）のラベルを持ったサブディレクトリがあるます。
その中にそれぞれの部位のDICOM画像が格納されています。（ラベル名に日本語使っても問題ありません）
ディレクトリの階層は次のようになります。

　　　　dcmdir2┬ Abdomen ┬ DICOM image a
　　　　　　　　│　　　　　│　　　・・・
　　　　　　　　│　　　　　└ DICOM image b
　　　　　　　　├ Chest　　┬ DICOM image c
　　　　　　　　│　　　　　│　　　・・・
　　　　　　　　│　　　　　└ DICOM image d
　　　　　　　　├ Head 　　┬ DICOM image e
　　　　　　　　│　　　　　│　　　・・・
　　　　　　　　│　　　　　└ DICOM image f
　　　　　　　　└ Lung  　　┬ DICOM image g
　　　　　　　　　　　　　　│　　　・・・
　　　　　　　　　　　　　　└ DICOM image h
　　　　　　　　　　　
実際にはタスクに応じてラベルは変更します。

プログラムの考え方

学習を行うとき，ラベル名は0からはじまるラベル番号に置き換えます。

プログラムのはじめに，ラベル名とラベル番号を対応させたファイルlabelName.csvを作成します。

画像をいくつのグループに分けるのか指定します。

例えば３つのグループ，0,1と2に分ける場合，この中から２つのグループを取り出して訓練データとしてTrainDataのDICOMファイル名のCSVファイルを作成します。
残りのひとつをTestDataとしてDICOMファイル名のCSVファイルを作成します。

グループに分割する方法は最も単純な方法を考えました。
３つに分割するのであれば，3で割ったあまりが，0か１か２でそれぞれのグループに分けます。

結果として３グループに分ける場合，訓練用とテスト用のファイルが次のように作成されます。
　　　　TrainData.0.3.csv　　　TestData0.3.csv
　　　　TrainData.１.3.csv　　　TestData１.3.csv
　　　　TrainData.2.3.csv　　　TestData2.3.csv
このCSVファイルには画像ファイル名，ラベル番号とラベル名がリストになっています。

プログラムを入力して実行してみよう

splitDicomData.py

import os

def saveLabelName(dataDir, saveDir):
    labelStr = ""
    labelName = {}
    dlist = os.listdir(dataDir)
    print("Subdirectory:",len(dlist),dlist)
    for i, subdir in enumerate(dlist):
        labelStr += str(i) + "," + subdir + "\n"
        labelName[subdir] = str(i)
   
    f = open(saveDir + "/labelName.csv", 'w')
    f.write(labelStr)
    f.close()
   
    return labelName

def splitData(dataDir, nGrp, saveDir, target, labelname):
    #リスト保存ファイル名
    saveTrainFname = saveDir + "TrainData." + str(target) + "." + str(nGrp) + ".csv"
    saveTestFname  = saveDir + "TestData."  + str(target) + "." + str(nGrp) + ".csv"
   
    strTrain = ""               #訓練画像のリストを保存する変数
    strTest  = ""               #テスト画像のリストを保存する変数
   
    dlist = os.listdir(dataDir) #カテゴリディレクトリの取得
    for subdir in dlist:
        flist = os.listdir(dataDir + "/" + subdir)
        print("\nLabel:",subdir,len(flist),"images")
        for j, fname in enumerate(flist):
            if j % nGrp == target:
                print(j,fname, " →　     Test")
                strTest += dataDir + "/" + subdir + "/"+fname +","
                strTest += labelName[subdir] + "," + subdir + "\n"
            else:
                print(j,fname, " →　Train")
                strTrain += dataDir + "/" + subdir + "/"+fname +","
                strTrain += labelName[subdir] + "," + subdir + "\n"
               
    #訓練データのリストファイルの保存           
    f = open(saveTrainFname,'w')
    f.write(strTrain)
    f.close()

    f = open(saveTestFname,'w')
    f.write(strTest)
    f.close()
   
if __name__=='__main__':
    dataDir = "../dcmdir2"           #対象画像の親ディレクトリ名
    nGroup = 3                       #グループ数
    saveDir = "./TrainTestList/"     #Listファイル保存ディレクトリ名
   
    labelName = saveLabelName(dataDir, saveDir)

    for target in range(nGroup):
        splitData(dataDir, nGroup, saveDir, target, labelName)

実行結果

プログラムを実行すると，はじめにラベル名とラベル番号のファイルが作成が行われます。
ラベル付けされたサブディレクトリの４つが表示されます。

$ python splitDicomData.py
Subdirectory: 4 ['Lung', 'Abdomen', 'Chest', 'Head']

結果として，labelName.csvがTrainTestListディレクトリ内に作成されます。
ファイルの内容は次のようにラベル番号とラベル名の対応表になっています。

続けてプログラムは，各サブディレクトリを読み，その中のDICOM画像を訓練用かテスト用か振り分けていきます。

Label: Lung 61 images
0 E5306265  →　     Test
1 W5305890  →　Train
2 K5303593  →　Train
3 M5305390  →　     Test
4 R5305640  →　Train
5 T5304000  →　Train
6 C5304421  →　     Test
7 A5304343  →　Train
8 D5304468  →　Train
9 O5303765  →　     Test
　　　　・・・・・
58 J5305250  →　Train
59 D5303265  →　Train
60 T5305734  →　     Test

Label: Abdomen 84 images
0 H5259375  →　     Test
1 F5258046  →　Train
2 L5258328  →　Train
　　　・・・・・

結果として，訓練用DICOMファイルのCSVファイルTrainData.i.n.csvとテスト用のTestData.i.n.csvが作成されます。
ファイルは，DICOM画像ファイル名，ラベル番号とラベル名のリストが保存されています。

このようにDICOM画像をディレクトリから移動させることなく，対象とする訓練データ，テストデータをファイルにリストとして保存することで柔軟にAIの実験を行うことができます。

プログラムの概説

saveLabelName関数

この関数はラベル名とラベル番号の対応をCSVファイルlabelName.csvを作成します。
２つの引数dataDirとsaveDirを持ちます。

• dataDirは，ラベル付けされたサブディレクトリの親ディレクトリです。
• saveDirはlabelName.csvをはじめ訓練・テスト用のDICOM画像リストのCSVファイルを格納するディレクトリ名です。

labelStrはCSVファイルにラベル番号とラベル名を書き出すための変数です。
labelNameはラベル名をキーにラベル番号を値にもつ辞書です。
これらを初期化しています。
　　　　labelStr = " "
　　　　labelName = { }

次にディレクトリdataDir内のディレクトリのリストをdlistに格納します。
　　　　dlist = os.listdir( dataDir )

次にfor文を使って，ディレクトリのリストdlistからひとつづつディレクトリ名をsubdirに取り出します。
それと0からの計数を変数iに入れ，ラベル番号（i）とラベル名（substr）のCSVデータを1行分作成しいます。
　　　　labelStr += str( i ) + " , " + subdir + " \n "
同時にラベル名（subdir）をキーにして，ラベル番号（i）を値に辞書lableNameに追加しています。
　　　　labelName[subdir] = str( i )

最後にCSVデータをファイルに書き出します。
　　　　f = open( saveDir + " /labelName.csv ",  'w ')
　　　　f.write( labelStr )
　　　　f.close( )
辞書データlabelNameを関数の返却値としてreturnしています。
この辞書は，次のsplitData関数で使います。

splitData関数

この関数は画像ファイルを訓練用，テスト用に分けます。
５つの引数を持ちます。

• dataDirは，DICOM画像が保存されている親ディレクトリ名です。
• nGrp,は，分割するグループ数です。
• saveDirは，訓練用・テスト用のDICOM画像リストファイルの保存ディレクトリ名です。
• targetは，テスト用の分割グループ番号です。指定された番号のグループがテスト用に，それ以外が訓練用にわけられます。
• labelnameは，ラベル名をキーにラベル番号を値に持つ辞書です。saveLabelName関数で作成されたものです。

ターゲットtargetとグループ番号nGrpから訓練用ファイル名saveTrainFnameとテスト用のファイル名saveTestFnameを作成しています。

訓練用とテスト用のCSVデータを作成するstrTrainとstrTestを初期化しています。

次に親ディレクトリdataDir内のサブディレクトリのリストをdlistに格納します。
for文によりひとつひとつディレクトリ名を取り出し，次の処理を行います。
取り出したひとつのサブディレクトリsubdir内のDICOM画像ファイルのリストflistを作成します。
　　　　flist = os.listdir( dataDir + " / " + subdir )

このファイルリストflistからfor文を用いて1画像づつ取り出し，その順序jも取り出します。
順序jが分割グールで割ったあまりが，もしターゲットtargetに等しければ，その画像をテスト用のCSVデータの変数に追加します。
　　　　if  j  %  nGrp  ==  target :
　　　　　　strTest += dataDir + " / " + subdir + " / "+fname +" , " 
さらに，ラベル番号を辞書labelNameを用いて取得して追加し，ラベル名subdirの値を追加します。

もし，順序jが分割グールで割ったあまりが，ターゲットtargetに等しけなければ，その画像を訓練用のCSVデータの変数に追加します。

最後に訓練用・テスト用のCSVデータをそれぞれのCSVファイルに保存して終了です。

main関数

この関数がプログラム実行時に最初に呼ばれ，

• 対象画像が入った親ディレクトリ
• 分割グループ数
• 各CSVファイルを保存するディレクトリ

が設定されています。

次にsaveLabelName関数を読んで，ラベル名の辞書をlabelNameに代入します。

グループ数分ループしてsplitData関数を呼び出して，訓練用・テスト用の画像リストのCSVファイルを作成します。

おわりに

このプログラムにより，簡単に訓練用・テスト用のDICOM画像のリストをCSVファイルに出力することができました。
また，このCSVファイルを編集することにより，容易に訓練用データやテスト用データを修正することができます。

次回，Pyゼミ3.04では今回作成したlabelName.csvや各訓練・テストCSVファイルを用いてCT画像の部位の分類を行ってみたいと思います。

課　題

今回，n枚の画像を分類するのに分割数mで割ったあまりを見て判断したが，乱数を用いてランダムにm個に分割するようにrandomSplitData関数を作成しなさい。
ただし，各グループの訓練用とテスト用にデータの重複があってはならない。

Pyゼミ3.02 手書き画像からDICOM画像へ

手書き数字画像MNISTから医療画像DICOMへ発展させよう

Keyword：MNIST，DICOM，Neural Network，Chainer

前回のPyゼミ3.01にて，手書き数字画像を用いて，ニューラルネットワークの学習を行ないました。
前回のニューラルネットワークをDICOM画像に発展させることを考えてみます。

MNISTがどのようなデータ構造をもった訓練画像なのか調べてみる

対話モードで，Chainerのdatasetsをimportして，訓練データ，テストデータをそれぞれtrainとtestに読み込んでみます。

>>> from chainer import datasets
>>> from chainer.datasets import tuple_dataset
>>>
>>> train, test = datasets.get_mnist(ndim = 3)

次に，訓練データの数をlenメソッドを使って調べてみます。

>>> len(train)
60000

確かに６万件のデータ（手書き数字画像）が読み込まれていることがわかります。
それでは，最初の０番めのデータについて，同じくいくつ格納されているのかlenメソッドで調べてみます。

>>> len(train[0])
2

２つです。
つまり，train[0][0]とtrain[0][1]の２つです。

さらにtrain[0][0]について，調べてみます。

>>> len(train[0][0])
1

１つです。
train[0][0][0]は１つということです。
さらに調べてみると，

>>> len(train[0][0][0])
28

28と表示されました。
これはMNISTの手書き数字の画素サイズです。
さらに調べると。

>>> len(train[0][0][0][0])
28

また，28と表示されました。


少しまとめてみましょう。
最初の60000と表示された。これはデータ数（画像と教師ラベル）を表しています。
つまり，train[0]からtrain[59999]のデータがあることを示しています。

先頭の0番目のデータについて注目してみると
この中には２つの情報があるようです。
ひとつめの情報（tarin[0][0]）を表示してみましょう。

>>>  train[0][0]
array([[[ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
          0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
　　　       ・・・略・・・
          0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
          0.        ,  0.        ,  0.        ]]], dtype=float32)

かなり大きなリストのデータです。多分画像じゃないかって想像がつきます。

ふたつめの情報(train[0][1])を表示してみましょう。

>>> train[0][1]

5

５と表示された。これはラベルデータじゃないかと想像できます。

次に，画像の方をさらに見てみると，len(train[0][0])は１だったので，train[0][0][0]しかない。

これを表示してみると，

>>> train[0][0][0]
array([[ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
　　　　　         ・・・略・・・
         0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ]], dtype=float32)

となる。arrayの中で [ [ で囲まれているので二次元のリストだ，つまり画像が格納されていることがわかります。

画像かどうか表示してみましょう。OpenCVで簡単に見ることができます。
（OpneCVのインストールはPyゼミ1.02参照）

>>> import cv2
>>>
>>> cv2.imshow('image', (train[0][0][0]))
>>> cv2.waitKey(0)

教師ラベルが５であったように，確かに５の手書き数字です。

少し補足するとtrain[0][0][0]は28×28画素のモノクロの画像1枚（１プレーン，１チャンネル）なので，len( train[0][0] ）は1を出力しました。
もし，カラー画像の場合，len( train[0][0] )は3を出力します。
つまり，RGBそれぞれがtrain[0][0][0]，train[0][0][1]，train[0][0][2]に格納されています。

まとめると

訓練画像データは5次元のデータで，ラベルデータは2次元である。
1次元：データの順序番号を表す。(0～データ数-1）
2次元：画像またはラベルを表す。0のとき画像，1のときラベル
3次元：画像プレーン番号（モノクロは0の１つのみ，カラーは0,1または2の３つ）
4次元：画像の行数（縦の画素数，画像の高さ）
5次元：画像の列数（横の画素数，画像の幅）

画像の画素x，yにアクセスするには次のようになります。
　　　　train[データ番号][0][プレーン番号][画像の行（y）][画像の列（ｘ）]
たとえば，10番めの画像のx＝11，y＝12の画素を得るには次のように書きます。
　　　　pix = train[9][0][0][12][11]

ラベルについては下記のようにアクセスできます。
　　　　tarin[データ番号][1]

したがって，画像とそのラベルが与えられたら5次元の訓練データを作成するプログラムを作成すれば良いことになります。

DICOM画像から訓練データを作成する

DICOM画像（このサンプルは個人の画像なのでこの実験以外では使用しないでください）から訓練データを作成する仕様としては以下のように考えます。

1. DICOM画像の入ったフォルダを指定して，その中の画像を訓練データとする。
2. ラベル名はファイル名（Brain01, Chest02など）から抽出し，それにラベル番号を自動で（0から）付与する（サンプル画像のファイル名の命名がラベル名となっている）。
3. DICOMファイル名，ラベル番号とラベル名のCSVファイルを作成する。
4. 上記，CSVファイルから訓練データを作成する関数を作る

プログラムの概要

constractTrainData関数

次の２つの引数を受け取とります。

• CSVファイル名（fname） ：DICOMファイル名，ラベル番号とラベル名のCSVファイル名を指定する。
• リサイズ値（resz）：画像を任意のサイズに拡大または縮小する画素数を指定する。

プログラムははじめに画像とラベルを格納するリストを初期化します。
　　　　image = [ ]
　　　　label = [ ] 

次にfor文にて，CSVファイルから1行づつ読み込みlineに代入し，読み込んだ回数を0から変数n代入します。
変数lineの値はカンマ区切りのデータで，data（リスト）に分割して格納します。
　　　　data[0]にはDICOMファイル名が
　　　　data[1]にはラベル番号が
　　　　data[2]にはラベル名が格納されます。

readDicom2png関数（2.03参照）を用いてウィンドウニングした画像をimgに代入します。
このとき，imgの画素値は0から255の8ビットの値をとります。

もし，リサイズの値reszが0より大きければリサイズしますが，0の場合はリサイズしません。

画像imgは浮動小数点32ビット（np.float32）に変換し，さらに0から255の値を-1.0から1.0の値に変換します。
　　　　img = (img - 128)/128

そして，imageのリストにappendメソッドを用いてimgを追加します。
　　　　image.append( [ img]  )

次にラベル番号は整数32ビット（int32）に変換して，labelのリストに追加します。
　　　　t = np.array( int( data[1] ), dtype = np.int32 )
　　　　label.append( t )

最後に，画像imageとラベルlabelのリストをTupleDatasetメソッドを使って，統合してtrainデータを作成します。
　　　　train = tuple_dataset.TupleDataset( image, label )

TupleDatasetは複数の同じ長さのデータセットをひとつのtuple型のデータセットにまとめるメソッドです。

プログラムを入力して実行してみよう

makeDicomTrain.py

import os, sys
import numpy as np
import cv2
import chainer
from chainer.datasets import tuple_dataset
import pydicom

def constractTrainData( fname , resz ):
    image = []
    label = []
    for n, line in enumerate(open(fname, 'r')):
        data = line.split(",")          # 画像ファイル名と正解データに分割
        print(n," File Name:", 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
        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      

if __name__=='__main__':

    dirname = "../dcmdir1/"
    fname   = "trainlist.csv"
    resize  = 200
    lblnum  = 0
    lblname = {}
    csvstr  = ""
   
    dlist = os.listdir(dirname)
    print(dlist)
    for nm in dlist:
        bsnm = nm[:-2]              #ファイル名の最後の2文字（数字）を取り除く
        if not bsnm in lblname:
            lblname[bsnm] = lblnum
            lblnum += 1
        csvstr += dirname+nm+","+str(lblname[bsnm])+","+bsnm+"\n"
    
    f = open(fname, 'w')
    f.write(csvstr)
    f.close()
      
    #画像データ，教師データの構築
    train = constractTrainData(fname, resize)
   
    print("\ntrain\n")
    print("画像データ数    :",len(train))
    print("画像/ラベル　　：",len(train[0]))
    print("画像チャンネル数:",len(train[0][0]))
    print("画像の高さ　　　:",len(train[0][0][0]))
    print("画像の幅　　　　:",len(train[0][0][0][0]))
    
    cv2.imshow("image",(train[0][0][0]))
    cv2.waitKey(0)
    print("\nLabel Number:",train[0][1])

実行結果

makeDicomTrain.pyを実行すると，指定したフォルダ（../dcmdir1/）内のファイルを表示します。
8枚のDICOM画像のファイル名とラベル番号が表示されます。
最後に訓練データの構成が表示されます。
画像をリサイズしたので画像の高さと幅が200となっています。

$ python makeDicomTrain.py
['Abdomen01', 'Chest02', 'Chest01', 'Brain01', 'Abdomen03', 'Abdomen02', 'Brain02', 'Chest03']

0  File Name: ../dcmdir1/Abdomen01  Label: 0
1  File Name: ../dcmdir1/Chest02  Label: 1
2  File Name: ../dcmdir1/Chest01  Label: 1
3  File Name: ../dcmdir1/Brain01  Label: 2
4  File Name: ../dcmdir1/Abdomen03  Label: 0
5  File Name: ../dcmdir1/Abdomen02  Label: 0
6  File Name: ../dcmdir1/Brain02  Label: 2
7  File Name: ../dcmdir1/Chest03  Label: 1

train

画像データ数    : 8
画像/ラベル　　： 2
画像チャンネル数: 1
画像の高さ　　　: 200
画像の幅　　　　: 200

Label Number : 0

trainlist.csvの出力

CSVファイルにDICOMファイル名，ラベル番号とラベル名が記録されてます。

確認画像のPNG出力と表示

tmp.pngファイルがreadDicom2png関数から出力されています。

mnistNN.pyにconstractTrainData関数を実装してみる

実装方法は以下の要領で行います。

1. mnistNN.pyをコピーしてdicomNN.pyに変更します。
2. makeDicomTrain.pyのconstractTrainData関数，readDicom2png関数をdicomNN.pyにコピーします。
3. mnistの入力部分は削除し，makeDicomTrain.pyのデータ入力部分を上記２つの関数のあとにコピーします。
4. os, sys, cv2とpydicomのインポートを追加します。
5. 今回，testデータはないので，trainデータの２つを使って実験をしてみましょう。

下記プログラムを参照してください。

dicomNN.py

import os, sys
import cv2
import numpy as np
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 n, line in enumerate(open(fname, 'r')):
        data = line.split(",")          # 画像ファイル名と正解データに分割
        print(n," File Name:", 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
        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   

#訓練データとテストデータの取得
dirname = "../dcmdir1/"
fname   = "trainlist.csv"
resize  = 28
lblnum  = 0
lblname = {}
csvstr  = ""

dlist = os.listdir(dirname)
print(dlist)
for nm in dlist:
    bsnm = nm[:-2]              #最後の2文字を取り除く
    if not bsnm in lblname:
        lblname[bsnm] = lblnum
        lblnum += 1
    csvstr += dirname+nm+","+str(lblname[bsnm])+","+bsnm+"\n"
f = open(fname, 'w')
f.write(csvstr)
f.close()

#画像データ，教師データの構築
train  = constractTrainData(fname, resize)
print("\ntrain\n")
print("画像データ数    :",len(train))
print("画像/ラベル　　：",len(train[0]))
print("画像チャンネル数:",len(train[0][0]))
print("画像の高さ　　　:",len(train[0][0][0]))
print("画像の幅　　　　:",len(train[0][0][0][0]))

class MyModel(Chain):                   #ネットワークモデルを定義する
    def __init__(self, nLabel):                 #ネットワークの定義
        super(MyModel, self).__init__(
            l1=L.Linear(784,100),       #全結合ネットワーク
            l2=L.Linear(100,100),           
            l3=L.Linear(100,nLabel),
        )
       
    def __call__(self, x):              #順伝播計算
        h1 = F.relu(self.l1(x))         #活性化関数にReLuを定義
        h2 = F.relu(self.l2(h1))
        return self.l3(h2)
       
model = MyModel(lblnum)                 #ネットワークモデルを生成する
model = L.Classifier(model, lossfun=F.softmax_cross_entropy)
optimizer = optimizers.Adam()           #最適化アルゴリズムを設定する
optimizer.setup(model)                  #モデルを最適化アルゴリズムにセットアップ

iterator = iterators.SerialIterator(train, 50, shuffle = True) #ミニバッチの設定
updater = training.StandardUpdater(iterator, optimizer) #更新処理の設定
trainer = training.Trainer(updater, (10, 'epoch'))      #訓練（エポック数）の設定
trainer.extend(extensions.ProgressBar())                #進捗表示の設定

trainer.run()                           #訓練実行

test=train[:2]  #試しに2つのデータについて

#訓練結果からtestデータを用いて推定を行う
ok = 0
n  = 0
for i in range(len(test)):              #テストデータ分繰り返す
    x = Variable(np.array([ test[i][0] ], dtype=np.float32)) #入力データ
    t = test[i][1]                      #正解データ
    y = model.predictor(x)              #テストデータから推定結果を計算
    out = F.softmax(y)
    ans = np.argmax(out.data)           #推定結果からラベルを得る
    n += 1
   
    np.set_printoptions(formatter={'float': '{: 0.2f}'.format})
    print('\n#',i)
    print('Infer :',out.data)           #推定結果を表示
    print('Sprvs :','  ','----->'*(int(t)),t)#教師データの表示
    if (ans == t):
        ok += 1   
        print('\tOK',ok/n)
     
print('平均正答率：', ok/len(test))

実行結果

mnistの手書きデータのサイズに合わせるために，DICOM画像を28×28にリサイズしています。
テストは未知のデータですべきですが，最初の２つの画像データでテストをしています。

今回はMNISTからDICOM画像への応用が成功しました。

$ python dicomNN.py
['Abdomen01', 'Chest02', 'Chest01', 'Brain01', 'Abdomen03', 'Abdomen02', 'Brain02', 'Chest03']
0  File Name: ../dcmdir1/Abdomen01  Label: 0
1  File Name: ../dcmdir1/Chest02  Label: 1
2  File Name: ../dcmdir1/Chest01  Label: 1
3  File Name: ../dcmdir1/Brain01  Label: 2
4  File Name: ../dcmdir1/Abdomen03  Label: 0
5  File Name: ../dcmdir1/Abdomen02  Label: 0
6  File Name: ../dcmdir1/Brain02  Label: 2
7  File Name: ../dcmdir1/Chest03  Label: 1

train

画像データ数    : 8
画像/ラベル　　： 2
画像チャンネル数: 1
画像の高さ　　　: 28
画像の幅　　　　: 28

# 0
Infer : [[ 0.67  0.22  0.11]]
Sprvs :     0
 OK 1.0

# 1
Infer : [[ 0.31  0.60  0.09]]
Sprvs :    -----> 1
 OK 1.0
平均正答率： 1.0

おわりに

MNISTはとても有名なデータセットなので，インターネットで調べると実はたくさん情報を得ることができます。
しかし，今回は解析的にデータ構造を追いかけることをしました。
こうした手法は，将来未知のデータセットに出会ったときに役に立つと思います。

今回訓練用データを作成するconstrauctTrain関数はこれからも使います。
また，画像ファイル名とラベル番号をCSVファイルに保存する考え方は，これからも訓練やテストのために使用します。

Pyゼミ2.02 PyQt5によるDICOM画像表示とページング

複数のDICOM画像をボタンをクリックしてページングしてみよう

keyword: PyQt5，DICOM，ページング

PyQt5のボタンの作り方とボタンをクリックした時の処理の仕方

前回（Pyゼミ2.01）はGUI上にDICOM画像を表示するだけでした。
今回は表示した画像の下に２つのボタン，NextボタンとPreviousボタンを付けます。

Nextボタンをクリックすると次のDICOM画像が表示され，
Previousボタンをクリックすると前のDICOM画像を表示します。

2つのボタンをQPushButtonクラスを使って次のように作成します。
        nextBtn = QPushButton("Next")          
        prevBtn = QPushButton("Previous")    

次にボタンをクリックした時(シグナルの発生）に行う処理（スロット）をclicked.connect関数を使って次のように記述します。
        nextBtn.clicked.connect(self.nextButtonClicked)            
        prevBtn.clicked.connect(self.prevButtonClicked)
nextBtnボタンがクリックされると，self.nextButtonClicked関数が呼びだされ，この関数に次の画像を読み込み，ウインドウに画像を表示するプログラムを記述します。

nextButtonClicked関数の処理

次の画像を表示するため画像ファイル番号fnoを一つインクリメントします。

        self.fno += 1

インクリメントしたfnoの値が画像枚数を超えたら，0にセットして元に戻す処理をしています。

        if self.fno == len(self.dcmfnm):

            self.fno = 0

DICOM画像のファイル名が保存されたリストdcmfnmのfno番目のDICOM画像を読み込み一旦PNGファイル（tmp.png）に保存します，

        self.readDicom2png(self.dirnm + self.dcmfnm[self.fno], 'tmp.png')

PNGファイルをQPixmapクラスでpixmapに読み込み，setPixmapクラスでラベルlblimgにセットすることで画像が更新されます。

        pixmap = QPixmap('tmp.png')

        self.lblimg.setPixmap(pixmap)
また，ファイル名も更新します。

        self.lblfnm.setText( "  File Name  : " + self.dcmfnm[ self.fno ] )

前回，GUIでDICOM画像を表示させるプログラムと比べて注意する点は以下のとおりです。

・ファイル名を表示するラベルlblfnmがself.lblfnmに変更しています。

・画像を表示するラベルlblimgがself.lblimgに変更しています。

これはself.を付けてインスタンス変数にすることで，クラス内で変数を共有することができます。
したがって，他の関数から値を変更したりすることができます。

プログラムを入力して実行してみましょう。

dicomPageButton.py

import sys,os
import numpy as np
import cv2
import pydicom 
from PyQt5.QtCore import *
from PyQt5.QtGui import *
from PyQt5.QtWidgets import *
from PyQt5.Qt import *

class DicomPageButton(QWidget):

    def __init__(self, dirname):
        super().__init__()
        self.dcmfnm = os.listdir(dirname)      # ディレクトリ内のDICOMファイル名
        self.dirnm  = dirname
        self.fno = 0    # 表示するファイル番号，最初は0番
        print("画像枚数：",len(self.dcmfnm))
        self.initUI()

    def initUI(self):
        #ファイル名を表示するラベルの作成
        self.lblfnm = QLabel("  File Name : " + self.dcmfnm[self.fno])
        fnmbox = QVBoxLayout()
        fnmbox.addWidget(self.lblfnm)

        #画像を読み込みラベルにセットする
        self.readDicom2png(self.dirnm + self.dcmfnm[self.fno], './tmp.png') 
        pixmap = QPixmap('./tmp.png')
        self.lblimg = QLabel()
        self.lblimg.setPixmap(pixmap)
        
        #画像表示レイアウトの設定
        imgbox = QVBoxLayout()
        imgbox.addWidget(self.lblimg)
        
        # NextとPreviousボタンの作成
        nextBtn = QPushButton("Next")          
        prevBtn = QPushButton("Previous")      
        nextBtn.clicked.connect(self.nextButtonClicked)            
        prevBtn.clicked.connect(self.prevButtonClicked)

        hbox = QHBoxLayout()       #ボタンを配置するための水平ボックス
        hbox.addWidget(prevBtn)
        hbox.addWidget(nextBtn)

        # 各パーツを配置する垂直なボックスを作成
        mainVbox = QVBoxLayout()        
        mainVbox.addLayout(fnmbox)
        mainVbox.addLayout(imgbox)
        mainVbox.addLayout(hbox) 

        # 垂直ボックスをウィンドウにレイアウトする
        self.setLayout(mainVbox)    
        self.setWindowTitle('DICOM Paging')    
        self.show()
        
    def nextButtonClicked(self):    #Nextボタンをクリックした時の処理
        self.fno += 1
        if self.fno == len(self.dcmfnm):
            self.fno = 0
        print("Push Next Button \t#", self.fno, "file:",self.dcmfnm[self.fno])
        self.readDicom2png(self.dirnm + self.dcmfnm[self.fno], 'tmp.png')
        self.pixmap = QPixmap('tmp.png')
        self.lblimg.setPixmap(self.pixmap)
        self.lblfnm.setText("  File Name  : " + self.dcmfnm[self.fno])


def prevButtonClicked(self):    #Priviousボタンをクリックした時の処理
        self.fno -= 1
        if self.fno == -1:
            self.fno = len(self.dcmfnm)-1
        print("Push Previous Button \t#", self.fno, "file:",self.dcmfnm[self.fno])
        self.readDicom2png(self.dirnm + self.dcmfnm[self.fno], 'tmp.png')
        self.pixmap = QPixmap('tmp.png')
        self.lblimg.setPixmap(self.pixmap)
        self.lblfnm.setText("  File Name : " + self.dcmfnm[self.fno])
                           
    #　DICOM画像読み込んでウィンドニング後PNGに保存       
    def readDicom2png(self, dcmfnm, tmpfnm):
        ds  = pydicom.read_file(dcmfnm)      #DICOM画像を読み込む
        wc  = ds.WindowCenter                #ウィンドウセンター値を代入
        ww  = ds.WindowWidth                 #ウィンドウ幅を代入        
        img = ds.pixel_array                 #画素値を代入               
        #表示画素値の最大と最小を計算する
        max = wc + ww / 2                      
        min = wc - ww / 2      
        #ウインドニング処理
        img = 255 * (img - min)/(max - min)   #最大と最小画素値を0から255に変換
        img[img > 255] = 255                  #255より大きい画素値は255に変換
        img[img < 0  ] = 0                    #0より小さい画素値は0に変換
        img = img.astype(np.uint8)
        cv2.imwrite(tmpfnm, img)

if __name__ == '__main__':
    app = QApplication(sys.argv)
    dirname = "../dcmdir1/" #DICOM画像の入ったディレクトリ
    ex = DicomPageButton(dirname)
    sys.exit(app.exec_())

実行結果

実行するとウィンドウが現れ，画像の下のボタンをクリックすると次の画像が表示されます。

ターミナルにはクリックスされたボタンの種類，ファイル番号とファイル名が表示されます。

$ python␣dicomPageButton.py⏎
画像枚数： 8
Push Next Button  # 1 file: Chest02
Push Next Button  # 2 file: Chest01
Push Next Button  # 3 file: Brain01

プログラムの説明

main関数

変数dirnameにDICOM画像が保存されたディレクトリを指定します。

これを引数にDicomPageButtonクラスのインスタンスを生成しています。

イニシャライザ__init__()関数

画像表示に必要な値の初期化を行っています。

DICOM画像の入ったディレクトリ（dirname）内のDICOMファイル名を読み込んでリストself.dcmfnmに代入しています。

        self.dcmfnm = os.listdir(dirname) 

DICOMファイル名も画像読み込み時に使用するので，self.dirnmに代入します。

        self.dirnm  = dirname

表示する画像番号self.fnoを初期化します。

initUI()関数

2つのラベルと2つのボタンを作成しています。

lblfnmは，DICOM画像のファイル名を表示するためのラベルです。

lblimgは，DICOM画像表示のためのラベルです。

nextBtnは，次の画像を表示するボタンです。

prevBtnは，前の画像を表示するボタンです。

各ラベルとボタンをレイアウトするBoxを作成しています。

fnmboxは，ラベルlblfnmをレイアウトします。

imgboxは，ラベルlblimgをレイアウトします。

hboxは，ボタンnextBtnとprevBtnを水平にレイアウトします。

mainVBoxは，上記2つのboxを垂直にレイアウトします。

最後にウィンドウに出力する処理を行います。

        self.setLayout(mainVbox)    

        self.setWindowTitle('DICOM Paging')    

        self.show()

nextButtonClicked()関数

self.dcmfnmの中の次のDICOM画像を読み込み，self.lblimgに表示する処理を行います。表示したDICOMファイル名を更新します．

prevButtonClicked()関数

self.dcmfnmの中のひとつ前のDICOM画像を読み込み，self.lblimgに表示する処理を行います。
表示したDICOMファイル名を更新します．

readDicom2png()関数

DICOM画像を読み込み，PNGファイルに出力します（Pyゼミ1.02参照）。

課　題

1）指定した画像番号のDICOMファイルとファイル名の更新の以下の4行が，iniUI関数の最初とnextButtonCliched関数とprevButtonClicked関数の3か所に記述されています．これをdisplayImageという関数にまとめなさい．
self.readDicom2png(self.dirnm + self.dcmfnm[self.fno], 'tmp.png') self.pixmap = QPixmap('tmp.png') self.lblimg.setPixmap(self.pixmap) self.lblfnm.setText(" File Name : " + self.dcmfnm[self.fno])

２）画像の表示順序がImage Positionに従っていません．頭部から腹部に向かってImage Positiondeソートして表示するように修正しなさい．（Pyゼミ1.04，dicomSort.pyを参照しなさい）

Pyゼミ1.02 DICOMからJPEGなどへの変換

DICOM画像をJPEGやPNGなどへの一般的な画像に変換しよう。

keyword: DICOM, OpenCV，JPEG, PNG，ImageJ

DICOM画像は専用のうDICOMビューワがないと表示することができません。
さいわいImageJのような高機能なフリーのビューワソフトも存在しますが，今回は一般的なJPEG画像など扱いやすい形に変換するプログラムを作成します。

プログラムはDICOM画像をウィンドニング処理（表示に最適な濃度変換処理）を加えてJPEGに変換します。
JPEG画像は非可逆な圧縮であり，圧縮時に画質を指定して保存することができます。
画質とファイルサイズの関係，画質低下に伴うエラーの状態を確認してみます。
さらに可逆圧縮（完全に元画像に復元できる圧縮方法）であるPNGファイルについても検討してみます。

OpenCVの画像ファイルの入出力関数を使う

OpneCV（Open Source Computer Vision Library）は画像処理をサポートするライブラリで，様々な関数が用意されて非常に便利なライブラリです。
このライブラリの中に画像の入出力もあるので，ここでは画像の書き込みにOpenCVのメソッドを使うことにします。
OpenCVの他のメソッドについては別途Pyゼミで紹介します。

OpenCVのインストールとバージョンの確認

OpenCVをインストールします。

$ pip␣install␣opencv-python⏎   

pythonを起動して，opneCVをインポートしてバージョンを確認してみます。

>>> import cv2
>>> cv2.__version__
'3.1.0'

OpenCVの3.1がインストールされているのが確認できました。

DICOM画像をJPEGに画像に変換して保存する

DICOM画像（この画像は実験にのみ使用してください）を読み込んでJPEGファイルに保存するプログラムを入力して実行してみましょう。

dicom2jpeg.py

import pydicom
import cv2
from matplotlib import pyplot as plt
    
def readDicom2jpeg(dcmfnm, savefnm):       #DICOM画像を読込，ウィンドニングしてJPEGに保存
    ds  = pydicom.read_file(dcmfnm)        #DICOM画像を読み込む
    wc  = ds.WindowCenter                　#ウィンドウセンター値を代入
    ww  = ds.WindowWidth                 　#ウィンドウ幅を代入       
    img = ds.pixel_array                 　#画素値を代入        
        
    #表示画素値の最大と最小を計算する
    max = wc + ww/2                     
    min = wc - ww/2      
    print("wc=",wc,"ww=",ww,"→ max=",max," min=",min)  
   
    #ウインドニング処理
    img = 255*(img - min)/(max - min)   　 #最大と最小画素値を0から255に変換        
    img[img > 255] = 255                  #255より大きい画素値は255に変換
    img[img < 0] = 0     　　　　　　　　 #JPEG画像として保存
    cv2.imwrite(savefnm, img, [cv2.IMWRITE_JPEG_QUALITY, 100])

if __name__ == '__main__':    
    readDicom2jpeg("../dcmdir1/Brain01","brain.jpg")   

このプログラムを実行するとカレントディレクトリにbrain.jpgが作成されます。

$ python␣dicom2jpeg.py⏎

生成されたbrain.jpgファイル

プログラムの解説

このプログラムは，プログラム内で使用するライブラリのインポート部とDICOM画像を読み込んでJPEG画像に保存するreadDicom2jpeg関数とプログラムの実行時に最初に呼び出されるmain関数から構成されています。

インポート部はDICOM画像を読み込むためのpydicomとJPEG画像を保存するためのOpenCVのインポートです。
後の実験のためにpyplotもインポートします。

readDicom2jpeg関数は，読み込むDICOMファイル名dcmfnmとJPEG画像を保存するファイル名savefnmを引数にとります。

はじめにDICOM画像を読み込み変数dsに代入します。

このdsからウィンドウセンタ，ウィンドウ幅と画像データをそれぞれ変数wc, ww, imgに代入しています。
imgにはCT画像の画素値が16ビットで格納されています。

次にウィンドニング処理を行います。
ウィンドニング処理はCT値の上限（max）を255に下限（min）を0になるように変換します。

maxとminはウィンドウセンターwcとウィンドウ幅wwから計算します。

    　　　max = wc + ww/2                   
    　　　min = wc - ww/2 　　　　

minとmaxの間のCT値を0から255の8bitの値に次式で変換します。

　　　　img = 255*(img - min)/(max - min)

従来（Javaなど）であれば，二重のfor文で各画素を変換するとろ，pythonは上のようにたった1行で処理してくれます。
Python素晴らしい。

変換後のimgの画素値が0未満の値（負の値）は0に，255より大きい値は255に変換します。

　　　　img[img > 255] = 255            

　　　　img[img < 0] =  0
この処理も従来（Javaなど）であれば，2重のfor文で各画素を変換していきますが，pythonでは上の2行で変換を行います。
しかも，for文で書くよりも高速に処理します。
Python素晴らしい。 

OpneCvのimwriteメソッドを使って，ウインドニング処理した画像imgをJPEG画像としてファイル名savefnmで保存します。

JPEG画像は劣化する

JPEG画像を保存するときに画質をcv2.IMWRITE_JPEG_QUALITYで指定することができます。
この値（Quality Factor, 以下QF）が100のとき，最高画質で保存します。
QFを小さくすると画質は低下しますが，圧縮したファイルサイズは小さくなります。

 それでは以下のプログラムをreadDicom2jpeg関数の最後に（JPEG保存のimwriteメソッドの下に）追加して実験してみましょう。
JPEGとは別に可逆圧縮（画質の劣化がない）のPNGでも保存してみましょう（ファイル名の拡張子を.pngにするだけ）。

    cv2.imwrite("QF70" + savefnm, img, [cv2.IMWRITE_JPEG_QUALITY, 70])
    cv2.imwrite("QF50" + savefnm, img, [cv2.IMWRITE_JPEG_QUALITY, 50])
    cv2.imwrite("QF30" + savefnm, img, [cv2.IMWRITE_JPEG_QUALITY, 30])
    cv2.imwrite("Brain.png", img)

作成されたJPEG画像のファイルサイズを見るとQFが小さくなると圧縮率も高くなり，ファイルサイズは小さくなっています。

　　　　　brain.jpg　　　   100.7kB
　　　　　QF70brain.jpg　　26.7kB
　　　　　QF50brain.jpg　　20.7kB
　　　　　QF30brain.jpg　　15.7kB

画像を拡大表示してみると

QFが70程度では画像の劣化は目立ちません（下図中）。
QFが30ではブロックノイズと呼ばれる8ｘ8画素の大きさのブロックが目立つようになります（下図右）。
JPEGは8ｘ8画素単位でDCT（discrete cosine transform，離散コサイン変換）を行い，高周波数成分を除くことで圧縮を行っています。
したがって，QFを小さくするとブッロク間の境界が目立つことになります（ブロックノイズ）。

オリジナルの画像との誤差を見てみよう

元のDICOMを読み込んだ画像データimgとJPEG圧縮した画像データimg2との差分により２つの画像の間の誤差を見てみましょう。

readDicom2jpeg関数の最後に下のプログラムを追加して実行してみよう。
imreadメソッドの引数の 0 は画像読み込みのフラグで，0はgray画像（1チャンネルの画像）として読み込むことを示します。
この引数はディフォルトでは3で，カラー画像として読み込みます。

    max = 5
    img2 = cv2.imread("QF30" + savefnm, 0)
    dif = abs(img - img2)
    plt.imshow(dif, cmap = 'gray', vmax = max, vmin = 0)
    plt.show()
    
    img2 = cv2.imread(savefnm, 0)
    dif = abs(img - img2)
    plt.imshow(dif, cmap = 'gray', vmax = max, vmin = 0)
    plt.show()
    
    img2 = cv2.imread("Brain.png", 0)
    dif = abs(img - img2)
    plt.imshow(dif, cmap = 'gray', vmax = max, vmin = 0)
    plt.show()

実行すると元のDICOM画像とQFを変えて保存した画像との差分画像（誤差画像）が表示されます。

QFが30との差分画像は明らかにブロックノイズが認められ，ノイズが多いのがわかります（下図左）。

QFを100にした最高品質のJPEG画像でも僅かにノイズがあるのがわかります（下図中）。

PNG画像は可逆圧縮であるため基本的に元画像との間に際はないはずだが，JPEGのQF100と同程度のノイズが認められます（下図右）。
これはOpenCVの内部的なアルゴリズムが原因かもしれません。

参考　ImageJによる差分画像

参考までにImageJで同様の実験結果を示します。
ImageJでDICOM画像を読み込み，8bitに変換した後にQFを変えたJPEG画像との差分（subtraction）を行いました。

QFが30のとき，差分画像は濃度変化の大きおエッジ部分に誤差が認められます。
QFが100のとき，わずかに誤差があるのが認められます。
PNG画像の差分画像には全く誤差が認められません。つまりPNGは完全に元に戻る画像圧縮フォーマットということです。

おわりに

今回はDICOM画像をウィンドニング処理を行い一般的な画像フォーマットに変換する方法を学びました。
この処理はこの後の様々なプログラムで再利用します。

Pyゼミ3.01 手書き認識のニューラルネットワーク

手書き認識のニューラルネットワークを試してみる

Keyword: MNIST, Neural Network, Chainer

だれもが通る道，MNISTの手書き数字の分類

誰もが簡単に取り組めるので，さまざまなサンプルが教科書やネットワーク上にあるので，既に経験したことがあるかもしれません。
手書き数字の分類のニューラルネットワークのプログラムを試しながら，医用画像の分類問題（Pyゼミ3.02）に進んでいきます。
もうすでに試した人には，復習です。

MNIST（Mixed National Institute of Standards and Technology database）は手書き数字画像を集めたデータベースで，訓練用に6万枚，テスト用に1万枚の手書き数字画像が用意されています。
このデータセットを使って誰もが簡単にニューラルネットワークの実験ができます。

プログラムの概要

ChainerではMNISTのデータを読み込みためのメソッドを用意しています。
　　　　train, test = datasets.get_mnist( ndim = 3 )
訓練用データとテスト用データがtrainとtestにそれぞれリスト形式で代入さています。

読み込まれたデータ数を確認するのにprint文で表示させてみると，
　　　　print( "train:", len( train ), "test:", len( test ) )
60000と10000と表示されます。

結構なデータ数なのでPCの性能によっては学習に時間がかかるかもしれません。
訓練用データを少し減らして実験するとよいです。
　　　　train = train[ : 600]　#先頭から600個の画像を訓練に使う
　　　　test  = test[ : 100]　　#先頭から100個の画像をテストに使う

ネットワークのモデルを作成する

MyModelクラス内でネットワークの定義を行います。
このネットワークは入力層（ノード数784）， ２つの隠れ層（ノード数100）と出力層（ノード数10）の単純なネットワークです。

init関数内（イニシャライザ）で各ネットワークを定義し，重みなどのパラメータの初期化を行います。

• l1は28×28画素（=784画素）の手書き数字画像を入力ノード数にもち，隠れ層（ノード数100）につなげるネットワークを定義しています。
• l2はノード数100の隠れ層から次のノード数100の隠れ層につなげるネットワークを定義しています。
• l3は前の隠れ層からの出力を10個のノード数（数字0から9の正解ラベルを学習する）を持つ出力層に接続するネットワークを定義しています。

call関数は順伝搬計算の手順を定義しています。

• 入力xをネットワークl1に与え，活性化関数ReLuを介して隠れ層の値h1を得る。
• 同様にh1をネットワークl2に与えて，h2を得る。
• 最後にh2をネットワークl3に与えた出力値（ベクトル）を返却する。

モデルの作成

先のネットワークの定義などMyModelクラスを呼び，モデルmodelを作成します。
　　　　model = MyModel( )
作成したモデルにおいて損失や精度を計算するためのClassifier関数を設定します。
このとき，損失関数の計算にSoftmax_Cross_Entropyを設定します。
　　　　model = L.Classifier( model, lossfun = F.softmax_cross_entropy )

次に，学習の最適化アルゴリズムAdamをoptimizerに設定して，作成したmodelにセットアップします。
　　　　optimizer = optimizers.Adam( )
　　　　optimizer.setup( model )
最適化アルゴリズムにはAdamの他にSGDなど様々なアルゴリズムが用意されています。

学習パラメータの設定と学習の実行

訓練データtrainを読み込むSerialIteratorに与えます。
引数にミニバッチサイズを50とデータセットのランダム化（shuffle）を有効にます。
shuffle = Trueにすることにより，毎エポックごとにシャッフルされてバッチデータが作成される。
　　　　iterator = iterators.SerialIterator( train, 50, shuffle = True )
iteratorに対して，optimaizerを設定して更新処理updaterを生成します。
　　　　updater = training.StandardUpdater( iterator, optimizer )
訓練のエポック数を設定してtrainerを生成する。ここでエポック数は10回です。
　　　　trainer = training.Trainer( updater, ( 10, 'epoch' ) )
学習の進捗をモニタするためにextendを使って進捗バーを表示するように設定します。
　　　　trainer.extend( extensions.ProgressBar( ) )
これらの一連の記述はお決まりです。
バッチサイズやエポック数を変えて実験することができます。

最後に生成したtrainerの実行を行います。
　　　　trainer.run()

このあと進捗バーが現れ，学習が行われます。
学習が終了すると，テスト用データを用いて推定を行います。
最終的に正答率を表示します。

プログラムを入力して実行してみよう

説明はこのくらいにしてプログラムを入力して実行してみよしょう。

mnistNN.py

import numpy as np
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

#訓練データとテストデータの取得
train, test = datasets.get_mnist(ndim=3) 
train = train[:6000]
test  = test[:1000]
print("train:", len(train), "test:", len(test))

class MyModel(Chain):                   #ネットワークモデルを定義する
    def __init__(self):                 #ネットワークの定義
        super(MyModel, self).__init__(
            l1=L.Linear(784,100),       #全結合ネットワーク
            l2=L.Linear(100,100),            
            l3=L.Linear(100,10),
        )
        
    def __call__(self, x):              #順伝播計算
        h1 = F.relu(self.l1(x))         #活性化関数にReLuを定義
        h2 = F.relu(self.l2(h1))
        return self.l3(h2)
        
model = MyModel()                       #ネットワークモデルを生成する
model = L.Classifier(model, lossfun=F.softmax_cross_entropy)
optimizer = optimizers.Adam()           #最適化アルゴリズムを設定する
optimizer.setup(model)                  #モデルを最適化アルゴリズムにセットアップ

iterator = iterators.SerialIterator(train, 50, shuffle = True) #ミニバッチの設定
updater = training.StandardUpdater(iterator, optimizer) #更新処理の設定
trainer = training.Trainer(updater, (10, 'epoch'))      #訓練（エポック数）の設定
trainer.extend(extensions.ProgressBar())                #進捗表示の設定

trainer.run()                           #訓練実行

#訓練結果からtestデータを用いて推定を行う
ok = 0
n  = 0
for i in range(len(test)):              #テストデータ分繰り返す
    x = Variable(np.array([ test[i][0] ], dtype=np.float32)) #入力データ
    t = test[i][1]                      #正解データ
    y = model.predictor(x)              #テストデータから推定結果を計算
    out = F.softmax(y)
    ans = np.argmax(out.data)           #推定結果からラベルを得る
    n += 1
    
    np.set_printoptions(formatter={'float': '{: 0.2f}'.format})
    print('\n#',i)
    print('Infer :',out.data)           #推定結果を表示
    print('Sprvs :','  ','----->'*(t),t)#教師データの表示
    if (ans == t):
        ok += 1    
        print('\tOK',ok/n)

        
print('平均正答率：', ok/len(test))

実行すると進捗バーが現れる。
エポック数，1秒間の繰り返し数や残り時間などの進捗状態が表示される。

$ python␣mnistNN.py⏎
train: 60000 test: 10000
     total [#################...............................] 37.50%
this epoch [###################################.............] 75.00%
      4500 iter, 3 epoch / 10 epoch
    406.63 iters/sec. Estimated time to finish: 0:00:18.444094

学習が終了すると，テスト用データで評価が表示されます。
各テスト用データ（#番号）について推定結果Inferのベクトルが表示されます。
この値が最も大きな値の要素がラベルの推定結果です。
その下の行に正解ラベルSprvsが----->のあとに表示されます。
一致していればOKの正答率は上がります。

    ・・・・・・・・・・・・

# 9995
Infer : [[ 0.00  0.00  1.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00]]
Sprvs :    ----->-----> 2
 OK 0.9770908363345339

# 9996
Infer : [[ 0.00  0.00  0.00  1.00  0.00  0.00  0.00  0.00  0.00  0.00]]
Sprvs :    ----->----->-----> 3
 OK 0.9770931279383815

# 9997
Infer : [[ 0.00  0.00  0.00  0.00  1.00  0.00  0.00  0.00  0.00  0.00]]
Sprvs :    ----->----->----->-----> 4
 OK 0.9770954190838168

# 9998
Infer : [[ 0.00  0.00  0.00  0.00  0.00  1.00  0.00  0.00  0.00  0.00]]
Sprvs :    ----->----->----->----->-----> 5
 OK 0.9770977097709771

# 9999
Infer : [[ 0.00  0.00  0.00  0.00  0.00  0.00  1.00  0.00  0.00  0.00]]
Sprvs :    ----->----->----->----->----->-----> 6
 OK 0.9771
正答率： 0.9771

最終的な正答率は0.9771でした。
非常に単純なネットワークであるが，良い結果が得られています。

プログラムの解説

インポート

Chainerを利用するのに必要なライブラリとnumpyをインポートします。

手書き数字データの準備

ChanerではMNISTのデータを読み込むメソッドが用意されています。
訓練データ６万件，テストデータ１万件がtrain,testそれぞれに格納されます。
引数のndim=3は数値画像の読み込み形式(1,28,28）の3次元を指定しています。
　　　　train, test = datasets.get_mnist(ndim=3) 

ネットワークモデルを定義するクラス

先に説明したように，init関数でネットワークを定義します。
call関数は入力xから順伝播計算を定義します。

モデルの生成と最適化アルゴリズムの設定

以下の4行はお決まりの文です。
　　　model = MyModel()                       #ネットワークモデルを生成する
　　　model = L.Classifier(model, lossfun = F.softmax_cross_entropy)
　　　optimizer = optimizers.Adam()     #最適化アルゴリズムを設定する
　　　optimizer.setup(model)      　　　#モデルを最適化アルゴリズムにセットアップ

学習パラメータの設定と実行

訓練データtrainを与え，バッチサイズ50で，シャッフルでバッチデータのランダム化を設定しています。
学習の最適化にoptimizerにAdamを設定しています。
実際の学習時にはこのStandardUpdaterにより重みの更新が行われます。
訓練回数を10エポックに設定しています。
学習の進捗をバー表示する設定を行い，最後にtrainerを実行します。


以上が学習のためのプログラムです。
学習が終了時には初期化されたl1,l2とl3の重みが更新され最適な重みになっています。
訓練結果の重みを利用して以下にtestのテストデータで正答率を求めてみましょう。

テストデータを用いいて手書き数字の推定を行う

推定に用いたテスト用データの数と推定結果が正解であった場合の数をそれぞれ変数nとokに代入するため初期化します。
次にfor文にてtestのデータ数分繰り返します。

入力ｘと正解tを準備し，推定結果を得る

まず，入力ｘと教師tを準備すします。
i番目の入力の手書き数字画像はtest[i][0]に格納されています。
これを入力ｘとしてmodelに与えるためにChainerのVariable関数を使って，32ビットの浮動小数点型のVariable変数に変換します。
　　　　x = Variable( np.array( [ test[ i ][ 0 ] ], dtype = np.float32 ) )
次に，i番目の正解データはtest[i][1]に格納されています。
これが正解tとなるため，そのまま代入します。
　　　　t = test[ i ][ 1 ] 
データの準備ができたので，modelのfwdInfer関数を使って入力データｘの推定を行います。
　　　　out = model.fwdInfer( x )
推定結果outはラベル０〜９に対する確率で，10次元のベクトルです。
このベクトルから最大値のラベルを求めて，推定結果ansを得ます。
　　　　ans = np.argmax( out.data )
変数ansには0～9までのいずれかの値が推定結果として格納されています。

推定結果と正答率の表示

浮動小数点表示を小数点以下２桁で表示するようにフォーマットの定義を行います。
　　　　np.set_printoptions(formatter={'float': '{: 0.2f}'.format})
データ番号を表示します。'\n'は改行を表す拡張表記です。
　　　　print( '\n#', i )
次に，推定結果のベクトルを表示します。
　　　　print( 'Infer :', out.data )          
その下に，正解ラベルを表示します。
正解ラベルtの数分，矢印を表示し，最後に正解ラベルtの値を表示します。
　　　　print( 'Sprvs :', '  ', '----->'*( t ), t )
もし，推定結果と正解ラベルが一致していれば，ベクトルの最大要素の下に正解ラベルが表示されるはずです。
最後に，if文で推定ラベルansと正解ラベルtを比較して一致したら変数okをインクリメントして，正答率を表示します。

課　題

今回，手書き数字画像のデータベースMNISTを使ってニューラルネットワークの学習と推定の実験を行った。
DICOM画像を読み込んでどのようにネットワークに与えて，学習と推定を行うか考えてみよう。

おわりに

MNISTはみな経験済みかもしれません。
これ以降，このプログラムをDICOM画像に発展させ，また畳み込みニューラルネットワークCNNに対応していきます。

Pyゼミ2.01 PyQt5によるDICOM画像表示．画像ファイル渡し

PyQt5を使ってDICOM画像を表示してみよう

Key word: PyQt5，　DICOM，

PyQtとは

PyQtは数あるGUI（Graphical User Interface）ライブラリ中でも人気があり，多くのサンプルプログラムなどがインターネット上に存在しています。
これらを使ってDICOM画像を表示するプログラムを作りながらPyQTの基本を学んでみようとおもいます。

PyQT5はAnacondaの中に含まれているのでAnacondaのインストール後にすぐに利用できます。

DICOM画像をPyQt5を使って表示する

DICOM画像を読み込み，いったんPNGファイルに保存した後，再びPyQTのQPixmapメソッドで読み込みウィンドウに表示します。
こうすることで，読込んだDICOMファイルを正しく読み込めたのかPNG画像で確認できます。
また，DICOM画像を使わないでPNGなどの一般的な画像への応用が容易になります。

DICOM画像はここにあります。実験以外では使用しないでください。

DicomDispayクラスを作成する

DicomDispayクラスは次の３つの関数で構成されている．
init関数：mainの中でDicomDisplayクラスを呼び出した際に初期化するイニシャライザです。
initUI関数：実際にGUIを作成する関数です。
これら２つの関数はPyQTでGUIを作成するための基本的な構成です。
そしてinitUI関数から様々なメソッドや関数を呼び出すことにより，目的に応じたGUIプログラムが作成することができます。

今回initUI関数がから呼び出される関数として，readDicom2png関数を作成します。
この関数はDICOMファイルを読み込みPNG画像に変換して保存する関数です。

プログラムの詳細は後にして，それではプログラムを入力して実行してみましょう。

DICOM画像表示プログラム

dicomDisplay.py

import sys,os
import numpy as np
import cv2
import pydicom
from PyQt5.QtCore import *
from PyQt5.QtGui import *
from PyQt5.QtWidgets import *
from PyQt5.Qt import *

class DicomDisplay(QWidget):

    def __init__(self, filename):
        super().__init__()
        self.dcmfnm = filename      # DICOMファイル名
        self.initUI()

    def initUI(self):
        #ファイル名を表示するラベルの作成
        lblfnm = QLabel(" File Name :  \t" + self.dcmfnm)
        fnmbox = QVBoxLayout()
        fnmbox.addWidget(lblfnm)

        #画像読み込みラベルにセットする
        self.readDicom2png(self.dcmfnm, './tmp.png')    #DICOMを読んでpngに出力
        pixmap = QPixmap('./tmp.png')
        lblimg = QLabel()
        lblimg.setPixmap(pixmap)
        print(">>pixmap type :", type(pixmap))
      
        #画像表示レイアウトの設定
        imgbox = QVBoxLayout()
        imgbox.addWidget(lblimg)

        # 各パーツを配置する垂直なボックスを作成
        mainVbox = QVBoxLayout()      
        mainVbox.addLayout(fnmbox)
        mainVbox.addLayout(imgbox)

        # 垂直ボックスをウィンドウにレイアウトする
        self.setLayout(mainVbox)  
        self.setWindowTitle('DICOM Display')  
        self.show()
              
      
    #　DICOM画像読み込んでウィンドニング後PNGに保存      
    def readDicom2png(self, dcmfnm, tmpfnm):
        ds  = pydicom.read_file(dcmfnm)      #DICOM画像を読み込む
        wc  = ds.WindowCenter                #ウィンドウセンター値を代入
        ww  = ds.WindowWidth                 #ウィンドウ幅を代入      
        img = ds.pixel_array                 #画素値を代入        
        #表示画素値の最大と最小を計算する
        max = wc + ww / 2                    
        min = wc - ww / 2    
        print("wc =",wc,"ww =",ww,"→ max =",max," min =",min)  
        #ウインドニング処理
        img = 255 * (img - min)/(max - min)   #最大と最小画素値を0から255に変換
        img[img > 255] = 255                  #255より大きい画素値は255に変換
        img[img < 0  ] = 0                    #0より小さい画素値は0に変換
        img = img.astype(np.uint8)
        print(">>pixel_array type:", img.dtype, type(img))
        cv2.imwrite(tmpfnm, img)

if __name__ == '__main__':
    app = QApplication(sys.argv)
    filename = "../dcmdir1/Brain01"
    ex = DicomDisplay(filename)
    sys.exit(app.exec_())

実行結果

実行するとウィンドウが現れ，DICOMファイル名とDICOM画像が表示されます。

$ python␣dicomDisplay.py⏎

プログラム解説

プログラムは実行すると最初に関連するライブラリをインポートしてmain関数から順に実行します。

main関数

最初に，QApplicationクラスはGUIを作成するための準備で最初に記述します。
filenameには表示するDICOM画像のファイル名を設定します。
次に，今回作成するDicomDisplayクラスをDICOMファイル名を引数に呼び出します。
最後のsys.exit( app.exec_( ) )はアプリケーションを終了するための処理を記述しています。

init関数は引数にselfとfilenameを受け取る。

この関数では，必要な変数などの初期化を行います。
表示するDICOM画像のファイル名を変数self.dcmfnmに代入します。
self.変数名はインスタンス変数と呼ばれ，クラス内で共有して利用できる便利な変数です。
実際にGraphical User Interfaceを作成するinitUI関数を呼び出します。

initUI関数は引数にselfを受け取ります。

この関数はselfを引数にしたので，イニシャライザで宣言したself.dcmfnmが関数内で利用可能です。

実際のGUIは２つのLabelと３つのVBoxLayoutを作成して以下の手順で行われます。

1. ファイル名を表示するためのLabel，lblfnmを作成する。
2. ラベルを配置するためのVBoxLayout，fnmboxを作成する。
3. 画像を表示するためのLabel，lblimgを作成する。
4. 画像のラベルを配置するためのVBoxLayout，imgboxの作成する。
5. ファイル名のfnmboxと画像のimgboxをを配置するVBoxLayoutとmainVboxを作成する。
6. 最後にmainVBoxをWindowにレイアウトして表示する。

という流れになります。

readDicom2png関数

この関数は引数にself, dcmfnmとtmpfnmの３つを受け取ります。
DICOMファイル（dcmfnm）を読み込み変数dsに代入します。
ウィンドニング処理をして最適な濃度に変換してPNGファイル（tmpfnm）に保存します。

ウィンドニング処理に必要なウィンドウセンターとウィンド幅の値をDICOMのタグから取得し変数wcとwwに代入します。
画素値のデータは変数imgに代入します。
画素値（ここではCT値）の最大値maxと最小値minをwcとwwから計算します。

実際のウィンドニング処理は次のように各画素に処理が行われます。
　　img = 255 * (img - min)/(max - min)
リストimgの値がを255を超えるとき次のように255に置き換えます。
　　img[img > 255] = 255
また，imgの値が負のとき，次のように0に置き換えます。
　　img[img < 0  ] = 0

最後にimgを符号なし8ビットに整えて，OpenCVのimwriteメソッドを使ってファイルに書き込みます。

おわりに

DICOM規格についてはPyゼミ1.01で説明しました。またDICOM画像のウィンドニング処理についても説明を行いましたので，ここでは説明が簡単になっています。

PyQt5は非常に簡単にGUIを作成することができます。
今後，ボタンの作成やマウスクリックのイベントの処理などを学びます。

DICOM画像を表示する際，いったんPNGのファイルに変換して保存してQPixmapに渡す，ファイル渡しの方法で画像表示しています。
filenameにPNGファイル名を指定して，QPixmapがこのファイルを読み込めば，PNG画像などにも容易に対応できることがわかると思います。

Pyゼミ1.01 DICOM画像表示

PythonでDICOM画像を表示する

キーワード：python, pydicom, pyplot, DICOM

医用画像の標準規格DICOMで保存された画像ファイルを読み込み表示してみる。

はじめてDICOM画像を表示したのはこのページを見てでした。
このプログラムがきっかけになり，Pythonに転向してしまいました。

Pydicomをインストールする

DICOM画像を読み込むにはPydicomのライブラリをインストールする必要があります。ターミナルから以下のようにしてpydicomをインストールします。
（␣は空白，⏎はenterキーです）

$pip␣install␣-U␣pydicom⏎

「-U」はインストール時にアップグレードを行うオプションです。
インストールが上手くいっているかpythonを起動して確認してみます。

>>>import pydicom
>>>from pydicom.data import get_testdata_files
>>>filename = get_testdata_files("rtplan.dcm")[0]
>>>ds = pydicom.dcmread(filename)  # plan dataset
>>>ds.PatientName
'Last^First^mid^pre'

Pydicomのサンプルを対話モードで入力してしてみます。
Pydicomインストール時にサンプルのDICOMファイル（”rtplan.dcm”)もインストールされているのでそのまま利用してみる。
（蛇足ですが，それにしてもrtplanとは放射線治療計画のDICOMデータなのだろうけど，ずいぶんマイナーなデータをサンプルに使ってますね）
変数dsにDICOMファイル”rtplan.dcm”の情報が格納されます。
ds.PatientNameはDICOMファイル内の患者氏名を指定している。ここでは匿名化された'Last^First^mid^pre'が表示されます。

デモ用DICOM画像のダウンロード

DICOM画像を読み込む準備をします。
デモ用のDICOM画像（私の過去のCT画像8枚です）をダウンロードし解凍します。
ディレクトリdcmdir1を現在のプログラムを作成するディレクトリ（例えばsrc01）の親ディレクトリ（一つ上のディレクトリ）にコピーします。

~/Documents/src01/readDicom.py 　#src01にプログラムが入る
~/Documents/dcmdir1/Brain01 　　#dcmdir1にDICOM画像が入る

DICOM画像を読み込んで表示する

それでは下記のプログラムをテキストエディタを使って入力してみます。
（テキストエディタにはgeditを使います。プログラムを作成するディレクトリ内で右クリックしリストから「新しいドキュメント/空のドキュメント」を指定し，「無題のドキュメント」を作成し，ファイルを右クリックして「名前の変更」を選択して，「readDicom.py」と変更します。）

readDicom.py

from matplotlib import pyplot as plt
import pydicom

d = pydicom.read_file('../dcmdir1/Brain01')

print(d)                #タグの表示
wc = d.WindowCenter     #ウィンドウセンターの取得
ww = d.WindowWidth      #ウィンドウ幅の取得
print("wc=", wc)
print("ww=", ww)

max = wc + ww / 2         #表示最大画素値
min = wc - ww / 2         #表示最小画素値
print("max=", max)
print("min=", min)

#画像の表示
plt.imshow(d.pixel_array, cmap = 'gray', vmax = max, vmin = min)
plt.show()

実行すると下のようにCT画像が表示されます。（ちょっと感動。はじめて表示したときは感動した。）

$python␣readDicom.py⏎

僅か10数行のプログラムでDICOM画像が表示さる。
こんなに簡単にDICOM画像が扱えることに感動。
私がPythonのとりこになった理由である。

プログラムの解説

はじめに，DICOM画像を処理できるようにpydicomライブラリをインポートし，画像表示のためのmatplotlibライブラリのpyplotモジュールをインポートする。

次にプログラムのあるディレクトリの親ディレクトリ（../）の下のdcmdir1の中のDICOMフィルBrain01を読み込み変数dに代入する。

変数dにはDICOMのメタ情報（患者氏名やID，検査日などの情報）と画像データが格納されている。変数dをprintするとメタ情報を表示できる。

(0008, 0000) Group Length                       UL: 430
(0008, 0008) Image Type                          CS: ['ORIGINAL', 'PRIMARY', 'AXIAL']
(0008, 0016) SOP Class UID                    UI: CT Image Storage
(0008, 0018) SOP Instance UID               UI: 1.2.392.200036.9116.2.6.1.16.161･････
(0008, 0020) Study Date                           DA: '20101112'
　　　　　　　　　　　　中略
(0010, 0010) Patient's Name                    PN: 'UESUGI^MASAHITO'
　　　　　　　　　　　　中略
(7005, 1030) [Main Modality in Study]         CS: 'CT'
(7005, 1040) [DLP Dose Length Product]   FD: 74.6
(7fe0, 0010) Pixel Data                　            OW: Array of 524288 bytes
wc= 40
ww= 80
max= 80.0
min= 0.0

DICOM画像を適正に表示するためにウィンドウセンターとウィンドウ幅をDICOMのメタ情報から取得wcとwwに代入する。
　　　　wc = d.WindowCenter     #ウィンドウセンターの取得
　　　　ww = d.WindowWidth      #ウィンドウ幅の取得

CT画像はCT値と呼ばれる値（画素値）で構成されていて，それは空気を-1000，水を0，骨を+1000として画像を構成している。したがって，CT装置のメーカーが異なってもCT値は同じように出力される。

そのまま表示するとコントラストの低い，メリハリのないボヤっとした画像になってしまう。
ウィンドウセンターとウィンドウ幅にはこの画像を表示するのに最適な条件が入っている。
ウィンドウセンターの値（デモ画像では40）を中心にウィンドウ幅の値（同80）で表示すると最適な表示条件で表示される。
最適な画像表示になるようにウィンドウセンターとウィンドウ幅を調整することをウィンドニングという。

最適な条件で表示するために，表示するCT値の上限maxと下限minを計算する。
上限max（プログラム結果では８０）よりも大きなCT値は白く，下限min（同０）よりも小さなCT値は黒く表示するようにする。

最大値maxはウィンドウ幅wwの半分の値をウィンドウセンターwcに加えて求める。
最大値maxより大きなCT値は白く，最小値minより小さなCT値は黒く表示する。
minとmaxで挟まれたCT値は0（黒）から255（白）の256諧調で表示されるため，コントラストのある見やすい画像が表示できる。

最後にimshowメソッド使って画像を表示する。
引数のd.pixel_arrayはDICOMファイルから読み込んだ画像データである。
cmap='gray'はモノクログレイ表示を指定している。
そしてvmax=max, vmin=minは画像データd.pixel_arrayの表示上限・下限範囲を指定した。
plt.show()メソッドで画像をウィンドウに表示する。
imshowメソッドの詳細はここから

まとめ

今回，pydicomライブラリを用いて簡単にDICOM画像を読込必要なメタ情報を表示した。さらに pyplotのメソッドを使ってDICOM画像を表示することができました。
これからpydicomを使いながら医用画像を人工知能学習させる方法（pyゼミ3で）を学びます。