Загрузка обучающих, валидационных и тестовых данных

В этом примере используется набор данных BraTS [2]. Набор данных BraTS содержит сканы МРТ опухолей головного мозга, а именно глиомы, которые являются наиболее распространенными первичными злокачественными новообразованиями головного мозга. Размер файла данных составляет ~ 7 ГБ. Если вы не хотите загружать набор данных BraTS, перейдите непосредственно в раздел Download Pretrained Network и Sample Test Set в этом примере.

Создайте директорию для хранения набора данных BraTS.

imageDir = fullfile(tempdir,'BraTS');
if ~exist(imageDir,'dir')
    mkdir(imageDir);
end

Чтобы скачать данные BraTS, перейдите на веб-сайт Medical Segmentation Decathlon и нажмите ссылку «Скачать данные». Загрузите файл «Task01_BrainTumour.tar» [3]. Разархивируйте файл TAR в директорию, заданную imageDir переменная. При успешном разархивировании imageDir будет содержать директорию с именем Task01_BrainTumour который имеет три подкаталога: imagesTr, imagesTs, и labelsTr.

Набор данных содержит 750 4-D томов, каждый из которых представляет собой стек 3-D изображений. У каждого 4-D объема есть размер 240 на 240 на 155 на 4, где первые три измерения соответствуют высоте, ширине и глубине 3D объемного изображения. Четвертая размерность соответствует различным модальностям скана. Набор данных разделен на 484 обучающих тома с вокселями и 266 тестовых томов, тестовые тома не имеют меток, поэтому этот пример не использует тестовые данные. Вместо этого пример разделяет 484 обучающих тома на три независимых набора, которые используются для обучения, валидации и проверки.

Предварительная обработка данных обучения и валидации

Чтобы обучить 3-D сеть U-Net более эффективно, предварительно обработайте данные МРТ с помощью вспомогательной функции preprocessBraTSdataset. Эта функция присоединена к примеру как вспомогательный файл.

Вспомогательная функция выполняет следующие операции:

Предварительная обработка данных может занять около 30 минут.

sourceDataLoc = [imageDir filesep 'Task01_BrainTumour'];
preprocessDataLoc = fullfile(tempdir,'BraTS','preprocessedDataset');
preprocessBraTSdataset(preprocessDataLoc,sourceDataLoc);

Создайте Datastore извлечения случайных закрашенных фигур для обучения и валидации

Используйте datastore случайного извлечения закрашенных фигур, чтобы передать обучающие данные в сеть и проверить процесс обучения. Этот datastore извлекает случайные закрашенные фигуры из основной истины изображений и соответствующих данных о пиксельных метках. Исправление является распространенным методом, чтобы предотвратить истощение памяти при обучении с произвольно большими томами.

Создайте imageDatastore для хранения 3-D данных изображения. Поскольку формат MAT-файла является нестандартным форматом изображения, необходимо использовать средство чтения MAT-файлов, чтобы включить чтение данных изображения. Можно использовать вспомогательную программу чтения MAT-файлов, matRead. Эта функция присоединена к примеру как вспомогательный файл.

volReader = @(x) matRead(x);
volLoc = fullfile(preprocessDataLoc,'imagesTr');
volds = imageDatastore(volLoc, ...
    'FileExtensions','.mat','ReadFcn',volReader);

Создайте pixelLabelDatastore для хранения меток.

lblLoc = fullfile(preprocessDataLoc,'labelsTr');
classNames = ["background","tumor"];
pixelLabelID = [0 1];
pxds = pixelLabelDatastore(lblLoc,classNames,pixelLabelID, ...
    'FileExtensions','.mat','ReadFcn',volReader);

Предварительный просмотр одного тома изображения и метки. Отображение маркированного объема с помощью labelvolshow функция. Сделать фон полностью прозрачным путем установки видимости метки фона (1) к 0.

volume = preview(volds);
label = preview(pxds);

viewPnl = uipanel(figure,'Title','Labeled Training Volume');
hPred = labelvolshow(label,volume(:,:,:,1),'Parent',viewPnl, ...
    'LabelColor',[0 0 0;1 0 0]);
hPred.LabelVisibility(1) = 0;

Создайте randomPatchExtractionDatastore который содержит обучающее изображение и данные о пиксельных метках. Определите размер закрашенной фигуры вокселей 132 на 132 на 132. Задайте 'PatchesPerImage' извлечь 16 случайным образом расположенных закрашенных фигур из каждой пары томов и меток во время обучения. Задайте мини-пакет размером 8.

patchSize = [132 132 132];
patchPerImage = 16;
miniBatchSize = 8;
patchds = randomPatchExtractionDatastore(volds,pxds,patchSize, ...
    'PatchesPerImage',patchPerImage);
patchds.MiniBatchSize = miniBatchSize;

Выполните те же шаги, чтобы создать randomPatchExtractionDatastore который содержит изображение валидации и данные о пиксельных метках. Можно использовать данные валидации, чтобы оценить, постоянно ли сеть учится, недофинансирована или избыточна по мере развития времени.

volLocVal = fullfile(preprocessDataLoc,'imagesVal');
voldsVal = imageDatastore(volLocVal, ...
    'FileExtensions','.mat','ReadFcn',volReader);

lblLocVal = fullfile(preprocessDataLoc,'labelsVal');
pxdsVal = pixelLabelDatastore(lblLocVal,classNames,pixelLabelID, ...
    'FileExtensions','.mat','ReadFcn',volReader);

dsVal = randomPatchExtractionDatastore(voldsVal,pxdsVal,patchSize, ...
    'PatchesPerImage',patchPerImage);
dsVal.MiniBatchSize = miniBatchSize;

Увеличение количества обучающих и валидационных данных при помощи transform функция с пользовательскими операциями предварительной обработки, заданными функцией helper augmentAndCrop3dPatch. Эта функция присоединена к примеру как вспомогательный файл.

The augmentAndCrop3dPatch функция выполняет следующие операции:

dataSource = 'Training';
dsTrain = transform(patchds,@(patchIn)augmentAndCrop3dPatch(patchIn,dataSource));

dataSource = 'Validation';
dsVal = transform(dsVal,@(patchIn)augmentAndCrop3dPatch(patchIn,dataSource));

Настройка 3-D слоев U-Net

В этом примере используется 3-D сеть U-Net [1]. В U-Net начальная серия сверточных слоев перемежается с максимальными слоями объединения, последовательно уменьшая разрешение входного изображения. Эти слои сопровождаются серией сверточных слоев, чередующихся с операторами повышающей дискретизации, последовательно увеличивая разрешение входного изображения. Перед каждым слоем ReLU вводят слой нормализации партии .. Имя U-Net происходит от того, что сеть может быть нарисована с симметричной формой, такой как буква U.

Создайте сеть 3-D U-Net по умолчанию с помощью unetLayers функция. Задайте сегментацию двух классов. Также задайте допустимое заполнение свертки, чтобы избежать программных продуктов при использовании стратегии перекрытия-плитки для предсказания тестовых объемов.

inputPatchSize = [132 132 132 4];
numClasses = 2;
[lgraph,outPatchSize] = unet3dLayers(inputPatchSize,numClasses,'ConvolutionPadding','valid');

Чтобы лучше сегментировать меньшие опухолевые области и уменьшить влияние больших фоновых областей, этот пример использует dicePixelClassificationLayer. Замените слой классификации пикселей слоем классификации пикселей Dice.

outputLayer = dicePixelClassificationLayer('Name','Output');
lgraph = replaceLayer(lgraph,'Segmentation-Layer',outputLayer);

Данные уже нормированы в разделе «Предварительная обработка данных обучения и валидации» этого примера. Нормализация данных в image3dInputLayer (Deep Learning Toolbox) не требуется, поэтому замените входной слой входным слоем, который не имеет нормализации данных.

inputLayer = image3dInputLayer(inputPatchSize,'Normalization','none','Name','ImageInputLayer');
lgraph = replaceLayer(lgraph,'ImageInputLayer',inputLayer);

Также можно изменить 3-D сеть U-Net с помощью Deep Network Designer App из Deep Learning Toolbox™.

Постройте график обновленной 3-D сети U-Net.

analyzeNetwork(lgraph)

Настройка опций обучения

Обучите сеть с помощью adam решатель оптимизации. Задайте установки гиперпараметров, используя trainingOptions (Deep Learning Toolbox) функция. Начальная скорость обучения устанавливается равным 5e-4 и постепенно уменьшается на протяжении всего периода обучения. Можно экспериментировать с MiniBatchSize свойство на основе памяти графический процессор. Чтобы максимизировать использование памяти графический процессор, предпочтите большие входные закрашенные фигуры для большого размера пакета. Обратите внимание, что слои нормализации партии . менее эффективны для меньших значений MiniBatchSize. Настройте начальную скорость обучения на основе MiniBatchSize.

options = trainingOptions('adam', ...
    'MaxEpochs',50, ...
    'InitialLearnRate',5e-4, ...
    'LearnRateSchedule','piecewise', ...
    'LearnRateDropPeriod',5, ...
    'LearnRateDropFactor',0.95, ...
    'ValidationData',dsVal, ...
    'ValidationFrequency',400, ...
    'Plots','training-progress', ...
    'Verbose',false, ...
    'MiniBatchSize',miniBatchSize);

Загрузите предварительно обученную сеть и выборку тестового набора

Загрузите предварительно обученную версию 3-D U-Net и пять выборочных тестовых томов и их соответствующих меток из набора данных BraTS [3]. Предварительно обученная данные модели и выборочные данные позволяют вам выполнить сегментацию тестовых данных, не загружая полный набор данных или ожидая обучения сети.

trained3DUnet_url = 'https://www.mathworks.com/supportfiles/vision/data/brainTumor3DUNetValid.mat';
sampleData_url = 'https://www.mathworks.com/supportfiles/vision/data/sampleBraTSTestSetValid.tar.gz';

imageDir = fullfile(tempdir,'BraTS');
if ~exist(imageDir,'dir')
    mkdir(imageDir);
end

downloadTrained3DUnetSampleData(trained3DUnet_url,sampleData_url,imageDir);

Обучите сеть

По умолчанию пример загружает предварительно обученную 3-D сеть U-Net. Предварительно обученная сеть позволяет запускать весь пример, не дожидаясь завершения обучения.

Чтобы обучить сеть, установите doTraining переменная в следующем коде, для true. Обучите модель с помощью trainNetwork (Deep Learning Toolbox) функция.

Обучите на графическом процессоре, если он доступен. Для использования GPU требуется Parallel Computing Toolbox™ и графический процессор с поддержкой CUDA ® NVIDIA ®. Для получения дополнительной информации смотрите Поддержку GPU by Release (Parallel Computing Toolbox). Обучение занимает около 30 часов на системе с несколькими графическими процессорами с 4 NVIDIA™ графическими процессорами Titan Xp и может занять еще больше времени в зависимости от вашего оборудования.

doTraining = false;
if doTraining
    modelDateTime = string(datetime('now','Format',"yyyy-MM-dd-HH-mm-ss"));
    [net,info] = trainNetwork(dsTrain,lgraph,options);
    save(strcat("trained3DUNet-",modelDateTime,"-Epoch-",num2str(options.MaxEpochs),".mat"),'net');
else
    inputPatchSize = [132 132 132 4];
    outPatchSize = [44 44 44 2];
    load(fullfile(imageDir,'trained3DUNet','brainTumor3DUNetValid.mat'));
end

Выполните сегментацию тестовых данных

Графический процессор настоятельно рекомендуется для выполнения семантической сегментации томов изображений (требует Parallel Computing Toolbox™).

Выберите источник тестовых данных, который содержит основную истину тома и метки для проверки. Если вы сохраняете useFullTestSet переменная в следующем коде как false, затем пример использует пять томов для проверки. Если вы задаете useFullTestSet переменная в trueзатем в примере используются 55 тестовых изображений, выбранных из полного набора данных.

useFullTestSet = false;
if useFullTestSet
    volLocTest = fullfile(preprocessDataLoc,'imagesTest');
    lblLocTest = fullfile(preprocessDataLoc,'labelsTest');
else
    volLocTest = fullfile(imageDir,'sampleBraTSTestSetValid','imagesTest');
    lblLocTest = fullfile(imageDir,'sampleBraTSTestSetValid','labelsTest');
    classNames = ["background","tumor"];
    pixelLabelID = [0 1];
end

The voldsTest переменная сохраняет основную истину тестовых изображений. The pxdsTest переменная сохраняет метки основной истины.

volReader = @(x) matRead(x);
voldsTest = imageDatastore(volLocTest, ...
    'FileExtensions','.mat','ReadFcn',volReader);
pxdsTest = pixelLabelDatastore(lblLocTest,classNames,pixelLabelID, ...
    'FileExtensions','.mat','ReadFcn',volReader);

Используйте стратегию перекрытия-плитки, чтобы предсказать метки для каждого тестового тома. Каждый тестовый том заполнен, чтобы сделать размер входа кратным выходу сигнала сети и компенсирует эффекты допустимой свертки. Алгоритм перекрытия-плитки выбирает перекрываемые закрашенные фигуры, предсказывает метки для каждой закрашенной фигуры при помощи semanticseg функция, а затем рекомбинирует закрашенные фигуры.

id = 1;
while hasdata(voldsTest)
    disp(['Processing test volume ' num2str(id)]);
    
    tempGroundTruth = read(pxdsTest);
    groundTruthLabels{id} = tempGroundTruth{1};
    vol{id} = read(voldsTest);
    
    % Use reflection padding for the test image. 
    % Avoid padding of different modalities.
    volSize = size(vol{id},(1:3));
    padSizePre  = (inputPatchSize(1:3)-outPatchSize(1:3))/2;
    padSizePost = (inputPatchSize(1:3)-outPatchSize(1:3))/2 + (outPatchSize(1:3)-mod(volSize,outPatchSize(1:3)));
    volPaddedPre = padarray(vol{id},padSizePre,'symmetric','pre');
    volPadded = padarray(volPaddedPre,padSizePost,'symmetric','post');
    [heightPad,widthPad,depthPad,~] = size(volPadded);
    [height,width,depth,~] = size(vol{id});
    
    tempSeg = categorical(zeros([height,width,depth],'uint8'),[0;1],classNames);
    
    % Overlap-tile strategy for segmentation of volumes.
    for k = 1:outPatchSize(3):depthPad-inputPatchSize(3)+1
        for j = 1:outPatchSize(2):widthPad-inputPatchSize(2)+1
            for i = 1:outPatchSize(1):heightPad-inputPatchSize(1)+1
                patch = volPadded( i:i+inputPatchSize(1)-1,...
                    j:j+inputPatchSize(2)-1,...
                    k:k+inputPatchSize(3)-1,:);
                patchSeg = semanticseg(patch,net);
                tempSeg(i:i+outPatchSize(1)-1, ...
                    j:j+outPatchSize(2)-1, ...
                    k:k+outPatchSize(3)-1) = patchSeg;
            end
        end
    end
    
    % Crop out the extra padded region.
    tempSeg = tempSeg(1:height,1:width,1:depth);

    % Save the predicted volume result.
    predictedLabels{id} = tempSeg;
    id=id+1;
end

Processing test volume 1
Processing test volume 2
Processing test volume 3
Processing test volume 4
Processing test volume 5

Сравнение основной истины с предсказанием сети

Выберите одно из тестовых изображений, чтобы оценить точность семантической сегментации. Извлеките первую модальность из объемных данных 4-D и сохраните этот объем 3-D в переменной vol3d.

volId = 1;
vol3d = vol{volId}(:,:,:,1);

Отобразите в формате montage центральный срез грунта и предсказанные метки в направлении глубины.

zID = size(vol3d,3)/2;
zSliceGT = labeloverlay(vol3d(:,:,zID),groundTruthLabels{volId}(:,:,zID));
zSlicePred = labeloverlay(vol3d(:,:,zID),predictedLabels{volId}(:,:,zID));

figure
montage({zSliceGT,zSlicePred},'Size',[1 2],'BorderSize',5) 
title('Labeled Ground Truth (Left) vs. Network Prediction (Right)')

Отображение помеченного основной истиной объема с помощью labelvolshow функция. Сделать фон полностью прозрачным путем установки видимости метки фона (1) к 0. Поскольку опухоль находится внутри ткани мозга, сделать некоторые воксели мозга прозрачными, чтобы опухоль была видна. Чтобы сделать некоторые воксели мозга прозрачными, задайте порог объема как число в области значений [0, 1]. Все нормированные интенсивности объема ниже этого порогового значения полностью прозрачны. Этот пример устанавливает порог объема меньше 1, чтобы некоторые пиксели мозга оставались видимыми, чтобы дать контекст пространственному расположению опухоли внутри мозга.

viewPnlTruth = uipanel(figure,'Title','Ground-Truth Labeled Volume');
hTruth = labelvolshow(groundTruthLabels{volId},vol3d,'Parent',viewPnlTruth, ...
    'LabelColor',[0 0 0;1 0 0],'VolumeThreshold',0.68);
hTruth.LabelVisibility(1) = 0;

Для того же объема отобразите предсказанные метки.

viewPnlPred = uipanel(figure,'Title','Predicted Labeled Volume');
hPred = labelvolshow(predictedLabels{volId},vol3d,'Parent',viewPnlPred, ...
    'LabelColor',[0 0 0;1 0 0],'VolumeThreshold',0.68);

hPred.LabelVisibility(1) = 0;

Это изображение показывает результат отображения срезов последовательно по всему объему. Помеченная основная истина находится слева, а сетевое предсказание - справа.

Количественная точность сегментации

Измерьте точность сегментации с помощью dice функция. Эта функция вычисляет коэффициент подобия Dice между предсказанной и основная истина сегментациями.

diceResult = zeros(length(voldsTest.Files),2);

for j = 1:length(vol)
    diceResult(j,:) = dice(groundTruthLabels{j},predictedLabels{j});
end

Вычислите средний счет Dice по набору тестовых томов.

meanDiceBackground = mean(diceResult(:,1));
disp(['Average Dice score of background across ',num2str(j), ...
    ' test volumes = ',num2str(meanDiceBackground)])

Average Dice score of background across 5 test volumes = 0.9993

meanDiceTumor = mean(diceResult(:,2));
disp(['Average Dice score of tumor across ',num2str(j), ...
    ' test volumes = ',num2str(meanDiceTumor)])

Average Dice score of tumor across 5 test volumes = 0.9585

Рисунок показывает boxplot (Statistics and Machine Learning Toolbox), который визуализирует статистику о счетах Dice по набору из пяти выборочных тестовых томов. Красные линии на графике показывают медианное значение Dice для классов. Верхняя и нижняя границы синей коробки обозначают 25-й и 75-й процентили соответственно. Черные усы распространяются на самые экстремальные точки данных, не считающиеся выбросами.

Если у вас есть Statistics and Machine Learning Toolbox™, то можно использовать boxplot функция для визуализации статистики о счетах Dice по всем вашим тестовым томам. Как создать boxplot, установите createBoxplot переменная в следующем коде, для true.

createBoxplot = false;
if createBoxplot
    figure
    boxplot(diceResult)
    title('Test Set Dice Accuracy')
    xticklabels(classNames)
    ylabel('Dice Coefficient')
end

Ссылки

[1] Чичек, В., А. Абдулкадир, С. С. Лиенкамп, Т. Брокс и О. Роннебергер. 3D U-Net: Обучение плотной объемной сегментации из разреженной аннотации ". В трудах Международной конференции по медицинскому Вычислению изображений и компьютерному вмешательству - MICCAI 2016. Афины, Греция, октябрь 2016, стр. 424-432.

[2] Isensee, F., P. Kickingereder, W. Wick, M. Bendszus и K. H. Maier-Hein. «Сегментация опухоли мозга и предсказание выживания радиомики: вклад в вызов BRATS 2017». В трудах BrainLes: Международный семинар по мозгу MICCAI. Квебек-Сити, Канада, сентябрь 2017, стр. 287-297.

[3] «Опухоли мозга». Медицинская сегментация десятичный атлон. http://medicaldecathlon.com/

Набор данных BraTS предоставляется Medical Segmentation Decathlon под лицензией CC-BY-SA 4.0. Все гарантии и заверения аннулируются; для получения дополнительной информации см. лицензию. MathWorks ® изменила набор данных, связанный в разделе Загрузка предварительно обученной сети и Набор тестов выборок этого примера. Модифицированный набор данных был обрезан до области, содержащей в основном мозг и опухоль, и каждый канал был нормализован независимо путем вычитания среднего значения и деления на стандартное отклонение обрезанной области мозга.

[4] Sudre, C. H., W. Li, T. Vercauteren, S. Ourselin, and M. J. Cardoso. Обобщенное перекрытие Dice как функция Глубокого обучения потерь для сильно несбалансированных сегментаций. Глубокое обучение в области анализа медицинских изображений и мультимодального обучения для поддержки клинических решений: Третий международный семинар. Квебек-Сити, Канада, сентябрь 2017, с. 240-248.

[5] Ronneberger, O., P. Fischer, and T. Brox. U-Net: Сверточные сети для сегментации биомедицинских изображений. В трудах Международной конференции по медицинскому Вычислению изображений и компьютерному вмешательству - MICCAI 2015. Мюнхен, Германия, октябрь 2015, с. 234-241. Доступно в arXiv:1505.04597.