文档

基于深度学习的三维脑肿瘤分割

这个例子展示了如何训练一个三维U-Net神经网络和执行语义分割脑肿瘤从三维医学图像。

该示例演示了如何训练一个三维U-Net网络,并提供了一个预训练网络。如果您选择训练网络,强烈建议使用支持cuda的NVIDIA™GPU,计算能力3.0或更高(需要并行计算工具箱™)。

介绍

语义分割包括用类标记图像中的每个像素或三维体素。

这个例子说明了在磁共振成像(MRI)扫描中使用深度学习方法对脑肿瘤进行语义分割。它还展示了如何执行二进制分割,其中每个体素都标记为肿瘤或背景。

医学图像分割的一个挑战是数据的类别不平衡,这妨碍了使用传统的交叉熵损失训练。这个例子通过使用加权多类骰子损失函数[4]来解决这个问题。对类进行权衡有助于抵消较大区域对Dice分数的影响,使网络更容易学习如何分割较小的区域。

本例使用3-D U-Net体系结构[1]执行脑肿瘤分割。U-Net是一种快速、高效且简单的网络,在语义分割领域非常流行。

下载培训、验证和测试数据

本例使用BraTS数据集[2]。BraTS数据集包含脑肿瘤的MRI扫描,即神经胶质瘤,这是最常见的原发性脑恶性肿瘤。数据文件大小为~ 7gb。如果您不想下载BraTS数据集,那么直接到下载预训练网络和样本测试集本例中的第节。

创建一个目录来存储BraTS数据集。

imageDir = fullfile (tempdir,“小鬼”);如果~存在(imageDir“dir”mkdir (imageDir);结束

要下载BraTS数据,请到医学分割Decathalon,然后点击“下载资料”链接。下载“task01_braintumor .tar”文件[3]。将TAR文件解压到指定的目录中imageDir变量。成功解压缩时,imageDir将包含一个名为Task01_BrainTumour它有三个子目录:imagesTrimagesTs标签

该数据集包含750个4-D卷,每个卷代表一堆3-D图像。每个4-D卷的尺寸为240 × 240 × 155 × 4,其中前三个维度对应于一个3-D体积图像的高度、宽度和深度。第四维对应不同的扫描方式。每个图像卷都有一个相应的像素标签。数据集分为484个训练卷和286个测试卷。

预处理培训和验证数据

为了更有效地训练三维U网络,使用辅助函数对MRI数据进行预处理preprocessBraTSdataset。此函数作为支持文件附加到示例中。万博1manbetx

helper函数执行以下操作:

  • 将数据裁剪到一个主要包含大脑和肿瘤的区域。裁剪数据减少了数据的大小,同时保留了每个MRI体积的最关键部分及其相应的标签。

  • 将每个体积的每个模态单独归一化,方法是减去平均值并除以剪裁后的大脑区域的标准偏差。

  • 将数据集分解为训练集、验证集和测试集。

数据预处理大约需要30分钟。

sourceDataLoc = [imageDir filesep .“任务01_脑瘤”];preprocessDataLoc = fullfile (tempdir,“小鬼”“preprocessedDataset”);preprocessBraTSdataset (preprocessDataLoc sourceDataLoc);

创建用于培训和验证的随机补丁提取数据存储

使用随机补丁提取数据存储将训练数据提供给网络,并验证训练进度。该数据存储从地面真实图像和相应的像素标签数据中提取随机补丁。补丁是一种常见的技术,用于在使用任意大容量训练时防止内存耗尽

首先,将训练图像存储在imageDatastore。由于MAT文件格式是非标准图像格式,因此必须使用MAT文件读取器才能读取图像数据。可以使用帮助器MAT文件读取器,matRead。此函数作为支持文件附加到示例中。万博1manbetx

volReader = @(x) matRead(x);volLoc = fullfile (preprocessDataLoc,“imagesTr”);volds = imageDatastore (volLoc,...“FileExtensions”“马特先生”“ReadFcn”, volReader);

创建一个pixelLabelDatastore储存标签。

labelReader=@(x)matRead(x);lblLoc=fullfile(预处理数据loc,“labelsTr”);一会= [“背景”“肿瘤”];pixelLabelID=[01];pixelLabelID pxds = pixelLabelDatastore (lblLoc,一会,...“FileExtensions”“马特先生”“ReadFcn”, labelReader);

预览一个图像卷和标签。显示标记卷使用拉伯沃尔秀函数。设置背景标签的可见性,使背景完全透明(1)到0

卷=预览(卷);标签=预览(pxds);图h=labelvolshow(标签,卷(:,:,:,1));h.LabelVisibility(1)=0;

创建一个randomPatchExtractionDatastore从图像数据存储和像素标签数据存储。指定64 × 64 × 64体素的patch大小。指定“PatchesPerImage”在训练过程中,从每对卷和标签中抽取16个随机定位的patch。指定迷你批处理大小为8。

patchSize=[64];patchPerImage=16;miniBatchSize=8;patchds=randomPatchExtractionDatastore(volds、pxds、patchSize、,...“PatchesPerImage”, patchPerImage);patchds。MiniBatchSize = MiniBatchSize;

通过使用使改变函数具有helper函数指定的自定义预处理操作增强3DPATCH.的增强3DPATCH函数随机旋转并反映训练数据,使训练更加稳健。此函数作为支持文件附加到示例中。万博1manbetx

dsTrain =变换(patchds @augment3dPatch);

创建一个randomPatchExtrationDatastore包含验证数据的。您可以使用验证数据评估网络是否随着时间的推移而不断学习、拟合不足或拟合过度。

volLocVal = fullfile (preprocessDataLoc,“imagesVal”);voldsVal = imageDatastore (volLocVal,...“FileExtensions”“马特先生”“ReadFcn”, volReader);lblLocVal = fullfile (preprocessDataLoc,“labelsVal”);pixelLabelID pxdsVal = pixelLabelDatastore (lblLocVal,一会,...“FileExtensions”“马特先生”“ReadFcn”, labelReader);dsVal = randomPatchExtractionDatastore (voldsVal pxdsVal patchSize,...“PatchesPerImage”, patchPerImage);dsVal。MiniBatchSize = MiniBatchSize;

设置3-D U-Net层

此示例使用三维U-Net网络的变体[1]。在U-Net中,最初的一系列卷积层穿插着最大池层,依次降低输入图像的分辨率。这些层之后是一系列卷积层,穿插着上采样运算符,依次提高输入图像的分辨率。U-Net的名称来源于fact表示网络可以用类似字母U的对称形状绘制。此示例修改U网络,在卷积中使用零填充,以便卷积的输入和输出具有相同的大小。

这个例子使用Deep Learning Toolbox™中的图层定义了3d U-Net,包括:

此示例还定义了一个自定义骰子丢失层,名为dicePixelClassification3dLayer,以解决数据[4]中的类不平衡问题。此层作为支持文件附加到示例中。有关详细信息,请参阅万博1manbetx定义自定义像素分类层与骰子损失

第一层,imageInput3dLayer,对大小为64×64×64体素的图像面片进行操作。

inputSize = [64 64 64 4];inputLayer = image3dInputLayer (inputSize,“归一化”“没有”“名字”“输入”);

图像输入层之后是三维U-Net的收缩路径。压缩路径由三个编码器模块组成。每个编码器包含两个带有3 × 3 × 3滤波器的卷积层,使特征图的数量增加一倍,然后使用reLu层进行非线性激活。第一个卷积之后是批处理归一化层。每个编码器以最大池化层结束,每个维度的图像分辨率减半。

给所有层赋予唯一的名称。编码器中的各层名称以“en”开头,后面跟着编码器模块的索引。例如,“en1”表示第一个编码器模块。

numFiltersEncoder = [32 64;64 128;128 256);层= [inputLayer];模块=1:3 modtag=num2str(模块);encoderModule=[Recolution3Dlayer(3,numFiltersEncoder(模块,1),...“填充”“一样”“增重剂”“窄正常”...“名字”,[“en”modtag,“_conv1”]);batchNormalizationLayer (“名字”,[“en”modtag,“_bn”]);reluLayer (“名字”,[“en”modtag,“_relu1”]);convolution3dLayer (numFiltersEncoder(模块,2),...“填充”“一样”“增重剂”“窄正常”...“名字”,[“en”modtag,“_conv2”]);reluLayer (“名字”,[“en”modtag,“_relu2”]);maxPooling3dLayer (2“步”2,“填充”“一样”...“名字”,[“en”modtag,“_maxpool”]);];层=[层;编码器模块];结束

创建三维U网的扩展路径。扩展路径由四个解码器模块组成。所有解码器都包含两个卷积层,带有3×3×3滤波器,可将特征映射的数量减半,然后使用reLu层进行非线性激活。前三个解码器以一个转置卷积层结束,该层对特征映射进行上采样最终解码器包括一个卷积层,将每个体素的特征向量映射到类。

为所有层指定唯一的名称。解码器中的层的名称以“de”开头,后跟解码器模块的索引。例如,“de4”表示第四个解码器模块。

numFiltersDecoder = [256 512;256 256;128 128;64 64);decoderModule4 =[卷积3dlayer (3,numFiltersDecoder(1,1),...“填充”“一样”“增重剂”“窄正常”...“名字”“de4_conv1”);雷卢耶(“名字”“de4_relu1”); 卷积3D层(3,numfilters解码器(1,2),...“填充”“一样”“增重剂”“窄正常”...“名字”“de4_conv2”);雷卢耶(“名字”“de4_relu2”);transposedConv3dLayer(2,numFiltersDecoder(1,2),“步”2,...“名字”“de4_transconv”);];decoderModule3=[卷积3Dlayer(3,numFiltersDecoder(2,1),...“填充”“一样”“增重剂”“窄正常”...“名字”“de3_conv1”);雷卢耶(“名字”“de3_relu1”);convolution3dLayer (numFiltersDecoder (2, 2),...“填充”“一样”“增重剂”“窄正常”...“名字”“de3_conv2”);雷卢耶(“名字”“de3_relu2”); transposedConv3dLayer(2,numFiltersDecoder(2,2),“步”2,...“名字”“de3_transconv”);];decoderModule2 =[卷积3dlayer (3,numFiltersDecoder(3,1),...“填充”“一样”“增重剂”“窄正常”...“名字”“de2_conv1”);雷卢耶(“名字”“de2_relu1”);卷积3层(3,numfilters解码器(3,2),...“填充”“一样”“增重剂”“窄正常”...“名字”“de2_conv2”);雷卢耶(“名字”“de2_relu2”);transposedConv3dLayer(2,numFiltersDecoder(3,2),“步”2,...“名字”“de2_transconv”);];

最后的解码器包括一个卷积层,该层将每个体素的特征向量映射到两个类(肿瘤和背景)。自定义Dice像素分类层权重损失函数,以增加小肿瘤区域对Dice得分的影响。

numLabels = 2;decoderModuleFinal =[卷积3dlayer (3,numFiltersDecoder(4,1),...“填充”“一样”“增重剂”“窄正常”...“名字”“de1_conv1”);雷卢耶(“名字”“de1_relu1”); 卷积3D层(3,numfilters解码器(4,2),...“填充”“一样”“增重剂”“窄正常”...“名字”“de1_conv2”);雷卢耶(“名字”“de1_relu2”);numLabels convolution3dLayer (1,“名字”“说服”);softmaxLayer(“名字”“softmax”);DicePixelsClassification3Dlayer(“输出”);];

将输入层和编码器模块与第四解码器模块连接。将其他解码器模块作为单独的分支添加到层图中。

层=[层;decoderModule4];lgraph = layerGraph(层);lgraph = addLayers (lgraph decoderModule3);lgraph = addLayers (lgraph decoderModule2);lgraph = addLayers (lgraph decoderModuleFinal);

使用级联层将每个编码器模块的第二reLu层与解码器模块中相同大小的转置卷积层连接起来。将每个级联层的输出连接到解码器模块的第一卷积层。

concat1 = concatenationLayer (4 2“名字”“concat1”);lgraph = addLayers (lgraph concat1);lgraph = connectLayers (lgraph,“en1_relu2”“concat1 /三机一体”);lgraph = connectLayers (lgraph,“de2_transconv”“concat1 / in2”);lgraph = connectLayers (lgraph,“concat1 /出”“de1_conv1”);concat2=连接层(4,2,“名字”“concat2”);lgraph = addLayers (lgraph concat2);lgraph = connectLayers (lgraph,“en2_relu2”“concat2/in1”);lgraph = connectLayers (lgraph,“de3_transconv”‘concat2/in2’);lgraph = connectLayers (lgraph,“concat2/out”“de2_conv1”);concat3=连接层(4,2,“名字”“concat3”);lgraph = addLayers (lgraph concat3);lgraph = connectLayers (lgraph,“en3_relu2”“concat3 /三机一体”);lgraph = connectLayers (lgraph,“de4_transconv”“concat3 / in2”);lgraph = connectLayers (lgraph,“concat3 /出”“de3_conv1”);

或者,您可以使用createUnet3dhelper功能创建三维U-Net层。这个函数作为支持文件附加到示例中。万博1manbetx

lgraph = createUnet3d (inputSize);

绘制图层图。

分析网络(lgraph)

指定培训选项

使用“亚当”优化求解器训练网络。属性指定超参数设置培训选项函数。初始学习率设置为5e-4,并随着训练的持续而逐渐降低。

选项=培训选项(“亚当”...“MaxEpochs”, 100,...“InitialLearnRate”5的军医,...“LearnRateSchedule”“分段”...“LearnRateDropPeriod”5,...“LearnRateDropFactor”, 0.95,...“验证数据”dsVal,...“验证频率”, 400,...“阴谋”“培训进度”...“详细”错误的...“MiniBatchSize”, miniBatchSize);

下载预训练网络和样本测试集

可选地,从BraTS数据集[3]下载预训练版本的3- d U-Net和5个样本测试卷及其相应的标签。预先训练的模型和样本数据使您能够对测试数据进行分割,而无需下载完整的数据集或等待网络训练。

trained3DUnet_url =“//www.tianjin-qmedu.com/万博1manbetxsupportfiles/vision/data/brainTumor3DUNet.mat”;sampleData_url ='//www.tianjin-qmedu.com/万博1manbetxsupportfiles/vision/data/sampleBraTSTestSet.tar.gz';imageDir = fullfile (tempdir,“小鬼”);如果~存在(imageDir“dir”mkdir (imageDir);结束downloadTrained3DUnetSampleData (trained3DUnet_url sampleData_url imageDir);

培训网络

配置培训选项和数据源后,使用trainNetwork作用要训练网络,请设置doTraining的参数符合事实的.具有CUDA功能的NVIDIA™ 强烈建议培训使用计算能力为3.0或更高的GPU。

如果你保留doTraining以下代码中的参数为,则该示例返回一个预先训练的3-D U-Net网络。

注:在NVIDIA上培训大约需要60小时™ 泰坦X和可能需要更长的时间取决于您的GPU硬件。

doTraining=false;如果doTraining model datetime=datestr(现在,“dd-mmm-yyyy-HH-MM-SS”);(网络,信息)= trainNetwork (dsTrain、lgraph选项);保存([“trained3DUNet -”modelDateTime的时代,num2str (maxEpochs)“马特先生”],“净”);其他的负载(fullfile (imageDir“trained3DUNet”“Braindunet.mat”));结束

你现在可以使用U-Net来对脑肿瘤进行语义分割。

执行测试数据的分段

选择包含地面真实体积和测试标签的测试数据源。如果保留useFullTestSet以下代码中的参数为,则示例使用5卷进行测试。如果你设置useFullTestSet参数到符合事实的,然后示例使用从完整数据集中选择的55个测试图像。

useFullTestSet = false;如果usfulltestset volLocTest = fullfile(preprocessDataLoc,“imagesTest”); lblLocTest=fullfile(预处理数据位置,“标签测试”);其他的volLocTest = fullfile (imageDir,“样本集”“imagesTest”);lblLocTest = fullfile (imageDir,“样本集”“标签测试”);一会= [“背景”“肿瘤”];pixelLabelID=[01];结束

使用helper函数将图像和标签的中心部分裁剪为尺寸为128 × 128 × 128的体素中央踏步机。此函数作为支持文件附加到示例中万博1manbetxvoldsTest变量存储地面真值测试图像。的PXD测试变量存储地面真相标签。

windowSize = [128 128 128];volReader = @(x) centerCropMatReader(x,windowSize);labelReader = @(x) centerCropMatReader(x,windowSize);voldsTest = imageDatastore (volLocTest,...“FileExtensions”“马特先生”“ReadFcn”,volReader);pXDTest=像素标签数据存储(lblLocTest、类名、像素标签ID、,...“FileExtensions”“马特先生”“ReadFcn”, labelReader);

对于每个测试图像,将地面真实图像的体积和标签添加到单元阵列。使用训练有素的网络semanticseg函数预测每个测试卷的标签。

在进行分割后,对预测标签进行后处理,将非脑体素标记为1,对应背景。使用测试图像来确定哪些体素不属于大脑。您还可以清除预测的标签,通过移除岛屿和填充洞使用medfilt3函数。medfilt3不支持分类数据,因万博1manbetx此将像素标签id转换为uint8然后,将过滤后的标签转换回分类数据类型,指定原始像素标签ID和类名。

id=1;虽然hasdata(卷测试)显示([“正在处理测试卷”num2str(id)])groundTruthLabels{id}=read(pxdsTest);vol{id}=read(voldsTest);tempSeg=semanticseg(vol{id},net);从测试图像中获取非大脑区域掩模。volMask=vol{id}(:,:,:,1)=0;%将预测标签的非脑区设置为背景。(1) tempSeg (volMask) =类名;%对预测的标签执行中值过滤。tempSeg = medfilt3 (uint8 (tempSeg) 1);%将筛选后的标签强制转换为“类别”。tempSeg=category(tempSeg,pixelLabelID,classNames);predictedLabels{id}=tempSeg;id=id+1;结束
加工测试第1卷加工测试第2卷加工测试第3卷加工测试第4卷加工测试第5卷

比较地面实况与网络预测

选择一幅测试图像来评估语义分割的准确性。从4-D体积数据中提取第一个模态,并将这个3-D体积存储在变量中vol3d

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

沿深度方向以蒙太奇方式显示地面真实和预测标签的中心切片。

zID =大小(vol3d, 3) / 2;zSliceGT = labeloverlay (vol3d (:,:, zID) groundTruthLabels {volId} (:,:, zID));zSlicePred = labeloverlay (vol3d (:,:, zID) predictedLabels {volId} (:,:, zID));图的标题(“标记地面真相(左)vs.网络预测(右)”)蒙太奇({zSliceGT; zSlicePred},“大小”(1 - 2),“BorderSize”5)

显示地面真相标记卷使用拉伯沃尔秀函数。设置背景标签的可见性,使背景完全透明(1)到0.因为肿瘤在脑组织内部,所以要使一些脑体素透明,这样肿瘤就可以被看到。为了使一些脑体素透明,将体积阈值指定为[0,1]范围内的一个数字。低于这个阈值的所有归一化体积强度都是完全透明的。本例将体积阈值设置为小于1,以便一些大脑像素保持可见,以提供肿瘤在大脑中的空间位置的上下文。

图h1 = labelvolshow(groundTruthLabels{volId},vol3d);h1.LabelVisibility (1) = 0;h1。VolumeThreshold = 0.68;

对于同一体积,显示预测的标签。

图h2 = labelvolshow(predictedLabels{volId},vol3d);h2.LabelVisibility (1) = 0;h2。VolumeThreshold = 0.68;

该图像显示了在整个卷中顺序显示切片的结果。

量化细分精度

利用骰子函数。该函数计算预测和地面真值分割之间的Dice相似系数。

diceResult = 0(长度(voldsTest.Files), 2);j=1:length(vol)diceResult(j,:)=dice(groundTruthLabels{j},predictedLabels{j});结束

计算一组测试卷的平均骰子分数。

meanDiceBackground =意味着(diceResult (: 1));disp ([“平均骰子得分的背景”,num2str(j),...' test volumes = 'num2str (meanDiceBackground)])
5个测试卷背景的平均骰子分数=0.99341
meanDiceTumor =意味着(diceResult (:, 2));disp ([“肿瘤的平均骰子得分”,num2str(j),...' test volumes = 'num2str (meanDiceTumor)])
5个测试体积中肿瘤的平均骰子分数=0.85204

图中显示了一个箱线图这可视化了关于五个样本测试卷集合中骰子分数的统计信息。图中的红线显示了类的骰子中值。蓝色框的上限和下限分别表示第25个和第75个百分位。黑色胡须延伸到未被视为异常值的最极端数据点。

如果您有统计学和机器学习工具箱™,那么您可以使用箱线图函数可视化关于所有测试卷的Dice分数的统计数据。创建一个箱线图,设置createBoxplot的参数符合事实的

createBoxplot = false;如果createBoxplot图形boxplot(dicerresult) title(“测试骰子精度”) xticklabels(类名)ylabel (“骰子系数”结束

总结

这个例子展示了如何创建和训练一个3-D U-Net网络来使用BraTS数据集执行3-D脑肿瘤分割。培训网络的步骤包括:

在训练三维U-Net网络或加载预训练的三维U-Net网络后,示例执行测试数据集的语义分割。该示例通过与地面真值分割进行视觉比较,并通过测量预测分割和地面真值分割之间的骰子相似系数来评估预测分割。

参考文献

[1] Cicek O。,A. Abdulkadir, S. S. Lienkamp, T. Brox, and O. Ronneberger. "3D U-Net: Learning Dense Volumetric Segmentation from Sparse Annotation." In国际医学图像计算和计算机辅助干预会议论文集希腊雅典,2016年10月,第424-432页。

[2] Isensee,F.,P.Kickingereder,W.Wick,M.Bendszus和K.H.Maier Hein.“脑肿瘤分割和放射组学生存预测:对BRATS 2017挑战的贡献”,摘自BrainLes会议录:国际MICCAI脑损伤研讨会.加拿大魁北克市,2017年9月,第287-297页。

[3]“脑癌”。医学分割十达隆。http://medicaldecathlon.com/

BraTS数据集由Medical Decathlon在4.0使用许可证。所有保证和声明均不受任何约束;有关详细信息,请参阅许可证。MathWorks®已修改了中链接的数据集下载预训练网络和样本测试集这个例子的一部分。修改后的样本数据集被裁剪到一个主要包含大脑和肿瘤的区域,每个通道通过减去平均值并除以裁剪后的大脑区域的标准差独立地进行归一化。

苏德雷,c.h., W. Li, T. Vercauteren, S. Ourselin, M. J. Cardoso。“广义骰子重叠作为高度不平衡分割的深度学习损失函数。”医学图像分析深度学习和临床决策支持的多模式学习:第三届国际研讨会万博1manbetx.加拿大魁北克市,2017年9月,第240-248页。

另请参阅

||||||

相关的例子

更多关于