MATLAB Examples

Contents

Medical Image Segmentation Using SegNet

本プログラムでは、MATLAB上でSegNetを構築・学習し、 学習済みネットワークを評価するところまでのワークフローを試行します。 画像データは血液塗抹標本画像を利用し、写っている寄生虫の部分、赤血球の部分と それ以外の領域の3クラスに分類することを目的とします。 画像データは米国CDC DPDx Parasite Image Libraryにて公開されているものを利用します。 https://www.cdc.gov/dpdx/malaria/index.html

This example shows how to create, train and evaluate SegNet for medical image segmentation. This example uses the images from CDC DPDx Parasite Image Linrary. https://www.cdc.gov/dpdx/malaria/index.html

clear all, close all, clc;

Setup

本デモではVGG-16をベースとしたSegNetを構築するため、 VGG-16を使えるように予めサポートパッケージがインストールされている必要があります。 インストール手順等につきましてはヘルプドキュメントをご覧ください。(>>doc('vgg16'))

vgg16();

SegNetのトレーニングにはGPUの利用を強く推奨します。 (GPUの利用にはParallel Computing Toolbox™が必要です。)

Load Blood Smear Images

血液塗抹標本画像をロードします。枚数はそれほど多くありませんが、 imageDatastoreを利用してワークスペースを効率的に利用します。

imgDir = fullfile(pwd,'images');
imds = imageDatastore(imgDir);

画像を1枚だけ読み込んで表示させます。

I = readimage(imds, 15);
imshow(I)

Label Ground Truth in a Collection of Images

血液塗抹標本画像に対し、ピクセルレベルでラベル付けを行います。 Image Labelerと呼ばれるアプリを利用することで、マウス操作でラベル画像の作成が行えます。

% groundTruthオブジェクトをワークスペースに作成
reconstLabelSession
labelDefs =

  3×3 table

          Name             Type       PixelLabelID
    ________________    __________    ____________

    'Background'        PixelLabel    [2]         
    'BloodParasites'    PixelLabel    [3]         
    'BloodCells'        PixelLabel    [1]         


labelData =

  16×1 table

                      PixelLabelData                  
    __________________________________________________

    'H:\BloodSmear\images\PixelLabelData\Label_01.png'
    'H:\BloodSmear\images\PixelLabelData\Label_02.png'
    'H:\BloodSmear\images\PixelLabelData\Label_03.png'
    'H:\BloodSmear\images\PixelLabelData\Label_04.png'
    'H:\BloodSmear\images\PixelLabelData\Label_05.png'
    'H:\BloodSmear\images\PixelLabelData\Label_06.png'
    'H:\BloodSmear\images\PixelLabelData\Label_07.png'
    'H:\BloodSmear\images\PixelLabelData\Label_08.png'
    'H:\BloodSmear\images\PixelLabelData\Label_09.png'
    'H:\BloodSmear\images\PixelLabelData\Label_10.png'
    'H:\BloodSmear\images\PixelLabelData\Label_11.png'
    'H:\BloodSmear\images\PixelLabelData\Label_12.png'
    'H:\BloodSmear\images\PixelLabelData\Label_13.png'
    'H:\BloodSmear\images\PixelLabelData\Label_14.png'
    'H:\BloodSmear\images\PixelLabelData\Label_15.png'
    'H:\BloodSmear\images\PixelLabelData\Label_16.png'


gTruth = 

  groundTruth with properties:

          DataSource: [1×1 groundTruthDataSource]
    LabelDefinitions: [3×3 table]
           LabelData: [16×1 table]

imageLabelerを起動し、予め作成してあるラベルを確認します。 imageLabeler起動 > Import Labels > From Workspace > gTruth選択 基本的に修正等を行う必要はなく、読み込んだラベル画像はそのまま以降の処理に 使うことができます。

sfile = fullfile(pwd, 'images', 'imageLabelingSession.mat');
imageLabeler(sfile)

Load Blood Smear Pixel-Labeled Images

今度はラベル画像を読み込みます。pixelLabelDatastoreを利用します。 塗抹標本のラベル画像は3個のクラスを持っていますが、ここでは2クラスに纏めて利用します。

classes = [
    "Background"
    "BloocCells"
    "BloodParasites"
    ];

pixelLabelID = cell(3,1);
pixelLabelID{1,1} = [2;0];
pixelLabelID{2,1} = 1;
pixelLabelID{3,1} = 3;

上で定義した2クラスとIDを利用し、pixelLabelDatastoreを作成します。

labelDir = fullfile(imgDir,'PixelLabelData');
pxds = pixelLabelDatastore(labelDir,classes,pixelLabelID);

ラベル画像を読み込み、該当する画像データにオーバーレイ表示させてみます。

C = readimage(pxds, 15);

cmap = bloodSmearColorMap;
B = labeloverlay(I,C,'ColorMap',cmap);

figure
imshow(B)
pixelLabelColorbar(cmap,classes);
% 何も色が付いていない領域はラベルを持っていない領域となり、
%  ネットワークのトレーニングには利用されません。

Analyze Dataset Statistics

各クラスラベルの分布を見るために、countEachLabelを利用します。 この関数を使うことで、クラス毎の総ピクセル数を確認することが出来ます。

tbl = countEachLabel(pxds)
tbl =

  3×3 table

          Name          PixelCount    ImagePixelCount
    ________________    __________    _______________

    'Background'        7.5845e+05    1.44e+06       
    'BloocCells'        6.0749e+05    1.44e+06       
    'BloodParasites'         74060    1.44e+06       

分布をFigureで可視化

frequency = tbl.PixelCount/sum(tbl.PixelCount);

figure
bar(1:numel(classes),frequency)
xticks(1:numel(classes))
xticklabels(tbl.Name)
xtickangle(45)
ylabel('Frequency')

理想的には、各クラスの総ピクセル数が等しくなっているのが望ましいですが、一般的には そのようなケースは稀となります。一部のクラスの総ピクセル数のみ極端に少ないような場合、 その領域が正しく検出できなくともAccuracyに大きく影響を与えないことになってしまいますので、 重み付けをしてから学習させます。

Resize Blood Smear Images

血液塗抹標本の画像サイズは300x300です。学習に要する時間やメモリ使用量を 削減したい場合はここでリサイズを実施します。

imageFolder = fullfile(imgDir,'imagesReszed',filesep);
imds = resizeBloodSmearImages(imds,imageFolder);

labelFolder = fullfile(imgDir,'labelsResized',filesep);
pxds = resizeBloodSmearPixelLabels(pxds,labelFolder);

Prepare Training and Test Sets

データセットの95%を学習に利用し、残りをテスト用とします。画像及びラベル画像を95:5に 分割します。

[imdsTrain, imdsTest, pxdsTrain, pxdsTest] = partitionBloodSmearData(imds,pxds);

分類した結果、学習用画像とテスト画像の総数がそれぞれどの程度あるか確認します。

numTrainingImages = numel(imdsTrain.Files)
numTestingImages = numel(imdsTest.Files)
numTrainingImages =

    15


numTestingImages =

     1

Create the Network

SegNetを作成します。Computer Vision System ToolboxにはsegnetLayersという関数があり 手軽にFCNを作成できるようになっています。 入力画像データのサイズ、クラス数を指定してSegNetを作成します。

imageSize = [300 300 3];
numClasses = numel(classes);
lgraph = segnetLayers(imageSize,numClasses,'vgg16');

Balance Classes Using Class Weighting

前述した通り各クラスの総ピクセル数が不均一であるため、総ピクセル数に応じた 重み付けを行います。

imageFreq = tbl.PixelCount ./ tbl.ImagePixelCount;
classWeights = median(imageFreq) ./ imageFreq
classWeights =

    0.8010
    1.0000
    8.2026

定義した重みをpixelClassificationLayerに反映させます。

pxLayer = pixelClassificationLayer('Name','labels','ClassNames', tbl.Name, 'ClassWeights', classWeights)
pxLayer = 

  PixelClassificationLayer with properties:

            Name: 'labels'
      ClassNames: {'Background'  'BloocCells'  'BloodParasites'}
    ClassWeights: [3×1 double]
      OutputSize: 'auto'

   Hyperparameters
    LossFunction: 'crossentropyex'

SegNetのpixelClassificationLayerを更新します。現在のpixelClassificationLayerを 削除し、新しく定義した層を追加します。

lgraph = removeLayers(lgraph, 'pixelLabels');
lgraph = addLayers(lgraph, pxLayer);
lgraph = connectLayers(lgraph, 'softmax' ,'labels');
%これで、SegNetの完成です。
figure, plot(lgraph)

Select Training Options

学習時のオプションを指定します。勾配法はMomentum-SGDを利用します。 利用するGPUに搭載されたメモリサイズに応じて、MiniBatchSizeを調整してください。

options = trainingOptions('sgdm', ...
    'Momentum', 0.9, ...
    'InitialLearnRate', 1e-3, ...
    'L2Regularization', 0.0005, ...
    'MaxEpochs', 3000, ...
    'MiniBatchSize', 1, ...
    'Shuffle', 'every-epoch', ...
    'Plots','training-progress', ...
    'VerboseFrequency', 1000);

Data Augmentation

学習に利用できる画像データの数が限られているため、Augmentation(数増し)を行います。 ここでは、ランダムに画像を回転させたり、X/Y方向に+/- 10ピクセルの範囲で平行移動させます。

augmenter = imageDataAugmenter('RandXReflection',true,'RandYReflection',true,...
    'RandXTranslation', [-10 10], 'RandYTranslation', [-10 10], 'RandRotation', [-180 180]);

Start Training

pixelLabelImageSourceを利用して、最終的に学習に用いるデータを定義します。 augmentationもここで行われます。

datasource = pixelLabelImageSource(imdsTrain,pxdsTrain,...
   'DataAugmentation',augmenter);

学習開始

[net, info] = trainNetwork(datasource,lgraph,options);
Training on single GPU.
Initializing image normalization.
|=========================================================================================|
|     Epoch    |   Iteration  | Time Elapsed |  Mini-batch  |  Mini-batch  | Base Learning|
|              |              |  (seconds)   |     Loss     |   Accuracy   |     Rate     |
|=========================================================================================|
|            1 |            1 |         1.29 |       1.2113 |       37.27% |       0.0010 |
|            3 |           45 |        29.06 |       1.1281 |       43.14% |       0.0010 |
|=========================================================================================|

Test Network on One Image

テスト用画像を1枚読み込み、結果を表示させます。トレーニングに利用した画像は 僅か十数枚でしたが、比較的綺麗に領域分割できていることが解ります。

idx = 1;
I = readimage(imdsTest,idx);
C = semanticseg(I, net);

% imshowで可視化
B = labeloverlay(I, C, 'Colormap', cmap, 'Transparency',0.4);
imshowpair(I, B, 'montage')
pixelLabelColorbar(cmap, classes);

セグメンテーションの結果と、真値(Ground truth)を比較してみます。 重ね書きした結果、緑やマゼンタになっている箇所が真値とは異なる箇所になります。

expectedResult = readimage(pxdsTest,idx);
actual = uint8(C);
expected = uint8(expectedResult);
imshowpair(actual, expected)

SegNetでセグメンテーションされた各クラスに対し、真値(Ground truth)領域が どの程度含まれているかを評価します。 これはIoUと呼ばれる指標で、jaccard関数を利用して測ることができます。

iou = jaccard(C, expectedResult);
table(classes,iou)
ans =

  3×2 table

        classes           iou  
    ________________    _______

    "Background"        0.44821
    "BloocCells"        0.13792
    "BloodParasites"    0.10171

真値に対する類似度を表現する係数としては、Jaccardの他にも Dice係数やBF係数などが良く用いられます。 MATLAB上ではそれぞれdice関数、bfscore関数を使って求めることができます。

Evaluate Trained Network

テスト用に切り分けておいた全ての画像を用い、SegNetの精度を測定します。

pxdsResults = semanticseg(imdsTest,net,'WriteLocation',tempdir,'Verbose',false);

semanticseg関数はpixelLabelDatastoreオブジェクトを返します。(ラベル画像のデータストア) 真値ラベルはpxdsTestデータストアにありますので、この2つのデータをevaluateSemanticSegmentation関数に与え、 各種評価指標を求めます。

metrics = evaluateSemanticSegmentation(pxdsResults,pxdsTest,'Verbose',false);

データセット全体に対する指標(平均精度、平均IoU値等)

metrics.DataSetMetrics
ans =

  1×5 table

    GlobalAccuracy    MeanAccuracy    MeanIoU    WeightedIoU    MeanBFScore
    ______________    ____________    _______    ___________    ___________

    0.46306           0.37906         0.22928    0.32682        0.2073     

各クラスに対する精度やIoU値

metrics.ClassMetrics
ans =

  3×3 table

                      Accuracy      IoU      MeanBFScore
                      ________    _______    ___________

    Background        0.60852     0.44821     0.32516   
    BloocCells        0.20127     0.13792     0.22616   
    BloodParasites    0.32739     0.10171    0.070578   

Supporting Functions

function pixelLabelColorbar(cmap, classNames)
%分類するクラスに対応するcolorbarを追加します。
% Add a colorbar to the current axis. The colorbar is formatted
% to display the class names with the color.

colormap(gca,cmap)

% Add colorbar to current figure.
c = colorbar('peer', gca);

% Use class names for tick marks.
c.TickLabels = classNames;
numClasses = size(cmap,1);

% Center tick labels.
c.Ticks = 1/(numClasses*2):1/numClasses:1;

% Remove tick mark.
c.TickLength = 0;
end
function cmap = bloodSmearColorMap()
% 血液塗抹標本画像の各クラスに対して紐付けられる色(colormap)を定義します。
% Define the colormap used by Blood Smear Pixel-Labeled Images.

cmap = [
    128 128 128   % Background
    000 000 192   % BloodCells
    128 000 000   % BloodParasites
    ];


% Normalize between [0 1].
cmap = cmap ./ 255;
end
function imds = resizeBloodSmearImages(imds, imageFolder)
% Resize images to [300 300].

if ~exist(imageFolder,'dir')
    mkdir(imageFolder)
else
    imds = imageDatastore(imageFolder);
    return; % Skip if images already resized
end

reset(imds)
while hasdata(imds)
    % Read an image.
    [I,info] = read(imds);

    % Resize image.
    I = imresize(I,[300 300]);

    % Write to disk.
    [~, filename, ext] = fileparts(info.Filename);
    imwrite(I,[imageFolder filename ext])
end

imds = imageDatastore(imageFolder);
end
function pxds = resizeBloodSmearPixelLabels(pxds, labelFolder)
% Resize pixel label data to [300 300].

classes = pxds.ClassNames;
labelIDs = 1:numel(classes);
if ~exist(labelFolder,'dir')
    mkdir(labelFolder)
else
    pxds = pixelLabelDatastore(labelFolder,classes,labelIDs);
    return; % Skip if images already resized
end

reset(pxds)
while hasdata(pxds)
    % Read the pixel data.
    [C,info] = read(pxds);

    % Convert from categorical to uint8.
    L = uint8(C);

    % Resize the data. Use 'nearest' interpolation to
    % preserve label IDs.
    L = imresize(L,[300 300],'nearest');

    % Write the data to disk.
    [~, filename, ext] = fileparts(info.Filename);
    imwrite(L,[labelFolder filename ext])
end

labelIDs = 1:numel(classes);
pxds = pixelLabelDatastore(labelFolder,classes,labelIDs);
end
function [imdsTrain, imdsTest, pxdsTrain, pxdsTest] = partitionBloodSmearData(imds,pxds)
% Partition Blood Smear Images by randomly selecting 95% of the data for training. The
% rest is used for testing.

% Set initial random state for example reproducibility.
rng(15);
numFiles = numel(imds.Files);
shuffledIndices = randperm(numFiles);

% Use 95% of the images for training.
N = round(0.95 * numFiles);
trainingIdx = shuffledIndices(1:N);

% Use the rest for testing.
testIdx = shuffledIndices(N+1:end);

% Create image datastores for training and test.
trainingImages = imds.Files(trainingIdx);
testImages = imds.Files(testIdx);
imdsTrain = imageDatastore(trainingImages);
imdsTest = imageDatastore(testImages);

% Extract class and label IDs info.
classes = pxds.ClassNames;
labelIDs = 1:numel(pxds.ClassNames);

% Create pixel label datastores for training and test.
trainingLabels = pxds.Files(trainingIdx);
testLabels = pxds.Files(testIdx);
pxdsTrain = pixelLabelDatastore(trainingLabels, classes, labelIDs);
pxdsTest = pixelLabelDatastore(testLabels, classes, labelIDs);
end