无监督医学图像去噪单元

这个例子使用了:

打开实时脚本

这个例子展示了如何使用UNIT神经网络从有噪声的低剂量CT图像生成高质量的计算机断层扫描(CT)图像。

注意:本示例引用了2021年5月1日访问的低剂量CT大挑战数据集。该示例使用来自数据集的胸部图像，这些图像现在处于限制访问之下。要运行此示例，您必须拥有与低剂量和高剂量CT图像兼容的数据集，并调整数据预处理和训练选项以适应您的数据。

这个例子使用了一个无监督的图像到图像转换(UNIT)神经网络，该神经网络从有限的数据样本中训练完整的图像。有关使用来自大量数据样本的图像数据块训练的CycleGAN神经网络的类似方法，请参见基于CycleGAN的无监督医学图像去噪(图像处理工具箱)．

x射线CT是一种流行的成像方式，用于临床和工业应用，因为它产生高质量的图像，并提供优越的诊断能力。为了保护病人的安全，临床医生建议使用低剂量的辐射。然而，低辐射剂量会导致图像的信噪比(SNR)降低，从而降低诊断的准确性。

深度学习技术为提高低剂量CT (LDCT)图像质量提供了解决方金宝搏官方网站案。使用生成对抗网络(GAN)进行图像到图像的转换，您可以将有噪声的LDCT图像转换为与常规剂量CT图像相同质量的图像。在此应用中，源域由LDCT图像组成，目标域由常规剂量图像组成。

CT图像去噪需要进行无监督训练的GAN，因为临床医生通常不会在同一时段获得同一患者的低剂量和常规剂量CT图像的匹配对。本例使用了支持无监督训练的UNIT体系结构。金宝app有关更多信息，请参见开始了解用于图像到图像转换的GANs(图像处理工具箱)．

下载LDCT数据集

本例使用来自低剂量CT大挑战的数据[2,3,4]。数据包括对常规剂量CT图像和模拟低剂量CT图像，包括99个头部扫描(标记N为神经)，100个胸部扫描(标记C为胸部)和100个腹部扫描(标记L为肝脏)。

指定dataDir作为数据集的期望位置。本例中的数据需要52 GB内存。

dataDir = fullfile(tempdir，“LDCT”，“LDCT-and-Projection-data”）;

要下载数据，请转到癌症影像档案的网站。这个例子只使用了两个病人的胸部扫描。从“Images (DICOM, 952 GB)”数据集中下载箱子文件“C081”和“C120”NBIA数据检索器．指定dataFolder变量作为下载数据的位置。下载成功后，dataFolder包含“C081”和“C120”两个子文件夹。

为培训、验证和测试创建数据存储

指定作为每个数据集的源的患者扫描。

scanDirTrain = fullfile(dataDir，“C120”，“08-30-2018-97899”）;scanDirTest = fullfile(dataDir，“C081”，“08-29-2018-10762”）;

创建imageDatastore管理用于训练和测试的低剂量和高剂量CT图像的对象。数据集由DICOM图像组成，因此使用自定义ReadFcn中的名称-值参数imageDatastore允许读取数据。

ext =“.dcm”；readFcn = @(x)rescale(dicomread(x));imdsLDTrain = imageDatastore(fullfile(scanDirTrain，"1.0000 -低剂量图像-71581"),.．.FileExtensions = ext ReadFcn = ReadFcn);imdsHDTrain = imageDatastore(fullfile(scanDirTrain，" 10000 -全剂量图像-34601"),.．.FileExtensions = ext ReadFcn = ReadFcn);imdsLDTest = imageDatastore(fullfile(scanDirTest，"1.000 -低剂量图像-32837"),.．.FileExtensions = ext ReadFcn = ReadFcn);imdsHDTest = imageDatastore(fullfile(scanDirTest，"1.0000 -全剂量图像-95670"),.．.FileExtensions = ext ReadFcn = ReadFcn);

从低剂量和高剂量CT训练数据集中预览训练图像。

lowDose =预览(imdsLDTrain);highDose = preview(imdsHDTrain);蒙太奇({lowDose, highDose})

预处理和增强训练数据

指定源和目标图像的图像输入大小。

inputSize = [256,256,1];

对训练数据进行增强和预处理变换属性指定的自定义预处理操作augmentDataForLD2HDCThelper函数。该函数作为支持文件附加到示例中。金宝app

的augmentDataForLD2HDCT函数执行以下操作:

使用双三次插值将图像调整为指定的输入大小。
在水平方向上随机翻转图像。
将图像缩放到范围[- 1,1]。这个范围与韵母的范围相匹配tanhLayer用于发电机。

imdsLDTrain = transform(imdsLDTrain， @(x)augmentDataForLD2HDCT(x,inputSize));imdsHDTrain = transform(imdsHDTrain， @(x)augmentDataForLD2HDCT(x,inputSize));

LDCT数据集提供了成对的低剂量和高剂量CT图像。然而，UNIT架构需要非配对数据进行无监督学习。本例通过在训练循环的每次迭代中变换数据来模拟未配对的训练和验证数据。

批训练和训练期间的验证数据

本例使用自定义训练循环。的minibatchqueue对象用于管理自定义训练循环中的观察数据的小批处理。的minibatchqueue对象也将数据强制转换为dlarray对象，用于在深度学习应用程序中实现自动区分。

指定小批量数据提取格式为SSCB(空间，空间，通道，批次)。设置DispatchInBackground返回的布尔值参数canUseGPU．如果支持的金宝appGPU可用于计算，则minibatchqueue对象在训练期间在并行池的后台预处理小批。

miniBatchSize = 1;mbqLDTrain = minibatchqueue(imdsLDTrain,MiniBatchSize= MiniBatchSize，.．.MiniBatchFormat =“SSCB”DispatchInBackground = canUseGPU);mbqHDTrain = minibatchqueue(imdsHDTrain,MiniBatchSize= MiniBatchSize，.．.MiniBatchFormat =“SSCB”DispatchInBackground = canUseGPU);

创建发电机网络

该单元由一个发生器和两个鉴别器组成。发生器从低剂量到高剂量进行图像到图像的转换。鉴别器是PatchGAN网络，它返回输入数据是真实的或生成的patch-wise概率。一个鉴别器区分真实和生成的低剂量图像，另一个鉴别器区分真实和生成的高剂量图像。

创建一个UNIT生成器网络unitGenerator(图像处理工具箱)函数。发生器的源编码器和目标编码器部分各由两个下采样块和五个剩余块组成。编码器部分共享五个剩余块中的两个。同样地，生成器的源和目标解码器部分各由两个下采样块和五个剩余块组成，解码器部分共享五个剩余块中的两个。

gen = unitGenerator(inputSize);

可视化发电机网络。

analyzeNetwork(创)

创建鉴别器网络

有两个鉴别器网络，每个图像域(低剂量CT和高剂量CT)一个。属性为源域和目标域创建标识符patchGANDiscriminator(图像处理工具箱)函数。

diskld = patchGANDiscriminator(inputSize,NumDownsamplingBlocks=4,FilterSize=3，.．.ConvolutionWeightsInitializer =“narrow-normal”NormalizationLayer =“没有”）;discHD = patchGANDiscriminator(inputSize，“NumDownsamplingBlocks”4 FilterSize = 3.．.ConvolutionWeightsInitializer =“narrow-normal”NormalizationLayer =“没有”）;

可视化鉴别器网络。

analyzeNetwork (discLD);analyzeNetwork (discHD);

定义模型梯度和损失函数

的modelGradientDisc而且modelGradientGen辅助函数分别计算鉴别器和生成器的梯度和损失。函数中定义了这些函数金宝app支持功能部分的示例。

每个鉴别器的目标是正确区分其域中图像的真实图像(1)和翻译图像(0)。每个鉴别器都有一个损失函数。

生成器的目标是生成经过翻译的图像，鉴别器将其分类为真实图像．发电机损耗是五种类型损耗的加权和:自重构损耗、循环一致性损耗、隐藏KL损耗、循环隐藏KL损耗和对抗损耗。

指定各种损失的权重因子。

lossWeights。自我ReconLossWeight = 10; lossWeights.hiddenKLLossWeight = 0.01; lossWeights.cycleConsisLossWeight = 10; lossWeights.cycleHiddenKLLossWeight = 0.01; lossWeights.advLossWeight = 1; lossWeights.discLossWeight = 0.5;

指定培训项目

指定Adam优化的选项。训练网络26个epoch。

numEpochs = 26;

为生成器和鉴别器网络指定相同的选项。

指定学习率为0.0001。
初始化尾随平均梯度和尾随平均梯度平方衰减率［］．
使用梯度衰减因子为0.5和平方梯度衰减因子为0.999。
使用系数为0.0001的权重衰减正则化。
使用1个小批量进行训练。

learnRate = 0.0001;gradDecay = 0.5;sqGradDecay = 0.999;weightDecay = 0.0001;genAvgGradient = [];genAvgGradientSq = [];discLDAvgGradient = [];discLDAvgGradientSq = [];discHDAvgGradient = [];discHDAvgGradientSq = [];

训练模型或下载预训练单元网络

默认情况下，该示例下载了NIH-AAPM-Mayo Clinic低剂量CT数据集的UNIT生成器的预训练版本。预训练的网络使您可以运行整个示例，而无需等待训练完成。

为了训练网络，设置doTraining变量转换为真正的．在自定义训练循环中训练模型。对于每个迭代:

方法读取当前小批处理的数据下一个函数。
方法计算鉴别器模型梯度dlfeval功能和modelGradientDischelper函数。
方法更新鉴别器网络的参数adamupdate函数。
方法评估生成器模型梯度dlfeval功能和modelGradientGenhelper函数。
方法更新发电机网络的参数adamupdate函数。
在每个纪元之后显示源域和目标域的输入和翻译图像。

如果有GPU，可以在GPU上进行训练。使用GPU需要并行计算工具箱™和支持CUDA®的NVIDIA®GPU。有关更多信息，请参见GPU计算要求(并行计算工具箱)．在NVIDIA Titan RTX上训练大约需要58个小时。

doTraining = false;如果doTraining创建一个显示结果的图形图(单位=“归一化”）;为iPlot = 1:4 ax(iPlot) = subplot(2,2,iPlot);结束迭代= 0;%遍历epoch为epoch = 1:numEpochs每个纪元洗牌数据重置(mbqLDTrain);洗牌(mbqLDTrain);重置(mbqHDTrain);洗牌(mbqHDTrain);运行循环，直到所有映像都在迷你批处理队列中% mbqLDTrain被处理而hasdata(mbqLDTrain)迭代=迭代+ 1;%从低剂量域读取数据imLowDose = next(mbqLDTrain);%从高剂量域读取数据如果hasdata(mbqHDTrain) == 0 reset(mbqHDTrain);洗牌(mbqHDTrain);结束imHighDose = next(mbqHDTrain);计算鉴别器梯度和损失[discldgrades, dischdgrades,discLDLoss,discHDLoss] = dlfeval(@modelGradientDisc，.．.创,discLD、discHD imLowDose、imHighDose lossWeights.discLossWeight);对低剂量鉴别器梯度应用权重衰减正则化discld= dlupdate(@(g,w) g+weightDecay*w, discld, discLD.Learnables);更新低剂量鉴别器参数[disclld,discLDAvgGradient,discLDAvgGradientSq] = adamupdate(disclld, discldgradients，.．.discLDAvgGradient discLDAvgGradientSq,迭代,learnRate、gradDecay sqGradDecay);对高剂量鉴别器梯度应用权重衰减正则化dischd= dlupdate(@(g,w) g+weightDecay*w, dischd, discHD.Learnables);%更新高剂量鉴别器参数[discHD,discHDAvgGradient,discHDAvgGradientSq] = adamupdate(discHD, dischdgradients, dischdgradientsq，.．.discHDAvgGradient discHDAvgGradientSq,迭代,learnRate、gradDecay sqGradDecay);计算发电机梯度和损耗[genGrad,genLoss,images] = dlfeval(@modelGradientGen,gen,discLD,discHD,imLowDose,imHighDose,lossWeights);对发电机梯度应用权重衰减正则化genGrad = dlupdate(@(g,w) g+weightDecay*w,genGrad,gen.Learnables);更新发电机参数[gen,genAvgGradient,genAvgGradientSq] = adamupdate(gen,genGrad,genAvgGradient，.．.迭代,genAvgGradientSq learnRate、gradDecay sqGradDecay);结束显示结果updateTrainingPlotLD2HDCT_UNIT (ax,图片{:});结束保存已训练的网络modelDateTime = string(datetime(“现在”格式=“yyyy-MM-dd-HH-mm-ss”));保存(fullfile (dataDir“trainedLowDoseHighDoseUNITGeneratorNet——”+ modelDateTime +“.mat”),“创”）;其他的net_url =“https://ssd.mathworks.com/金宝appsupportfiles/vision/data/trainedLowDoseHighDoseUNITGeneratorNet.zip”；downloadTrainedNetwork (net_url dataDir);负载(fullfile (dataDir“trainedLowDoseHighDoseUNITGeneratorNet.mat”));结束

利用训练网络生成高剂量图像

从低剂量测试图像的数据存储中读取并显示图像。

idxToTest = 1;imLowDoseTest = readimage(imdsLDTest,idxToTest);图imshow (imLowDoseTest)

将图像转换为数据类型单．将图像数据重新缩放到生成器网络的最后一层所期望的范围[- 1,1]。

imLowDoseTest = im2single(imLowDoseTest);imLowDoseTestRescaled = (imLowDoseTest-0.5)/0.5;

创建一个dlarray对象，该对象向生成器输入数据。如果支持的金宝appGPU可用于计算，则通过将数据转换为a在GPU上执行推理gpuArray对象。

dlLowDoseImage = dlarray(imlowdosetestrescale，“SSCB”）;如果canUseGPU dlLowDoseImage = gpuArray(dlLowDoseImage);结束

函数将输入的低剂量图像转换为高剂量域unitPredict(图像处理工具箱)函数。生成的图像的像素值范围为[- 1,1]。为了显示，将激活缩放到范围[0,1]。

dlImLowDoseToHighDose = unitPredict(gen,dlLowDoseImage);imHighDoseGenerated = extractdata(gather(dlImLowDoseToHighDose));imHighDoseGenerated = rescale(imHighDoseGenerated);imshow (imHighDoseGenerated)

读取并显示地面真实高剂量图像。高剂量和低剂量测试数据存储不被打乱，因此地面真实高剂量图像直接对应于低剂量测试图像。

imHighDoseGroundTruth = readimage(imdsHDTest,idxToTest);imshow (imHighDoseGroundTruth)

以蒙太奇的方式显示输入的低剂量CT、生成的高剂量版本和真实的高剂量图像。尽管该网络是在单个患者扫描数据上进行训练的，但该网络可以很好地测试来自其他患者扫描的图像。

imshow([imLowDoseTest imHighDoseGenerated imHighDoseGroundTruth]) title(“测试图像”+ num2str (idxToTest) +"低剂量，产生的高剂量，和真实的高剂量"）

金宝app支持功能

模型梯度函数

的modelGradientGenHelper函数计算生成器的梯度和损失。

函数[genGrad,genLoss,images] = modelGradientGen(gen,discLD,discHD,imLD,imHD,lossWeights) [imLD2LD,imHD2LD,imLD2HD,imHD2HD] = forward(gen,imLD,imHD);hidden = forward(gen,imLD,imHD,Outputs=“encoderSharedBlock”）;[~，imLD2HD2LD,imHD2LD2HD，~] = forward(gen,imHD2LD,imLD2HD);cycle_hidden = forward(gen,imHD2LD,imLD2HD,Outputs=“encoderSharedBlock”）;计算不同损失selfReconLoss = computeReconLoss(imLD,imLD2LD) + computeReconLoss(imHD,imHD2HD);hiddenKLLoss = computeKLLoss(隐藏);cycleReconLoss = computeReconLoss(imLD,imLD2HD2LD) + computeReconLoss(imHD,imHD2LD2HD);cycleHiddenKLLoss = computeKLLoss(cycle_hidden);outA = forward(discLD,imHD2LD);outB = forward(discHD,imLD2HD);advLoss = computeAdvLoss(outA) + computeAdvLoss(outB);将发电机的总损耗计算为五个损耗的加权和genTotalLoss =.．.selfReconLoss * lossWeights。selfReconLossWeight +.．.hiddenKLLoss * lossWeights。hiddenKLLossWeight +.．.cycleReconLoss * lossWeights。cycleConsisLossWeight +.．.cycleHiddenKLLoss * lossWeights。cycleHiddenKLLossWeight +.．.advLoss * lossWeights.advLossWeight;更新发电机参数genGrad = dlgradient(genTotalLoss,gen.Learnables);将数据类型从单数组转换为单数组genLoss = extractdata(genTotalLoss);images = {imLD,imLD2HD,imHD,imHD2LD};结束

的modelGradientDiscHelper函数计算两个鉴别器的梯度和损失。

函数[discldgrades, dischdgrades,discLDLoss,discHDLoss] = modelGradientDisc(gen，.．.discLD,discHD,imRealLD,imRealHD,discLossWeight) [~，imFakeLD,imFakeHD，~] = forward(gen,imRealLD,imRealHD);计算低剂量图像鉴别器的损耗outRealLD = forward(discLD,imRealLD);outFakeLD = forward(discLD,imFakeLD);discLDLoss = disclosight *computeDiscLoss(outRealLD,outFakeLD);更新低剂量图像鉴别器参数discldgrades = dlgradient(discLDLoss,discLD.Learnables);计算高剂量图像鉴别器的损耗outRealHD = forward(discHD,imRealHD);outFakeHD = forward(discHD,imFakeHD);discHDLoss = disclosight *computeDiscLoss(outRealHD,outFakeHD);%更新高剂量图像鉴别器参数dischdgrades = dlgradient(discHDLoss,discHD.Learnables);将数据类型从单数组转换为单数组discLDLoss = extractdata(discLDLoss);discHDLoss = extractdata(discHDLoss);结束

损失函数

的computeDiscLossHelper函数计算鉴别器损失。每个鉴别器损失是两个分量的和:

1向量与鉴别器在真实图像上的预测值之间的平方差， $Y_{真正的}$
零向量与鉴别器对生成图像的预测之间的平方差， ${\hat{Y}}_{翻译}$

$discriminatorLoss ＝ {（ 1 - Y_{真正的} ）}^{2} + {（ 0 - {\hat{Y}}_{翻译} ）}^{2}$

函数discLoss = computeDiscLoss(Yreal,Ytranslated) discLoss = mean(((1-Yreal).^2)，“所有”) +.．.意思是(((0-Ytranslated) ^ 2),“所有”）;结束

的computeAdvLoss辅助函数计算发电机的对抗损失。对抗损失是1的向量与翻译图像上的鉴别器预测之间的平方差。

$adversarialLoss ＝ {（ 1 - {\hat{Y}}_{翻译} ）}^{2}$

函数advLoss = computeAdvLoss(Ytranslated) advLoss = mean(((Ytranslated-1).^2)，“所有”）;结束

的computeReconLoss辅助函数计算发生器的自重构损失和周期一致性损失。自我重建的损失是 $l^{1}$ 输入图像与自重建图像之间的距离。循环一致性损失是 $l^{1}$ 输入图像与其循环重建版本之间的距离。

$selfReconstructionLoss ＝ {为（ Y_{真正的} - Y_{自我 - 重建} ）为}_{1}$

$cycleConsistencyLoss ＝ {为（ Y_{真正的} - Y_{周期 - 重建} ）为}_{1}$

函数reconLoss = computeReconLoss(Yreal,Yrecon) = mean(abs(Yreal-Yrecon)，“所有”）;结束

的computeKLLosshelper函数计算发电机的隐藏KL损失和循环隐藏KL损失。隐藏KL损失是零向量和'encoderSharedBlock'对自我重建流的激活。周期隐藏KL损失是零向量和'encoderSharedBlock'循环重建流的激活。

$hiddenKLLoss ＝ {（ 0 - Y_{encoderSharedBlockActivation} ）}^{2}$

$cycleHiddenKLLoss ＝ {（ 0 - Y_{encoderSharedBlockActivation} ）}^{2}$

函数klLoss = computeKLLoss(hidden) klLoss = mean(abs(hidden.^2)，“所有”）;结束

参考文献

[1]刘明宇，Thomas Breuel, Jan Kautz，“无监督图像到图像的翻译网络”。在神经信息处理系统研究进展，2017.https://arxiv.org/pdf/1703.00848.pdf．

[2] McCollough, C.H.， Chen, B.， Holmes, D.， III .，段晓霞，Yu, Z.， Yu, L.，冷，S.， Fletcher, J.(2020)。数据来自低剂量CT图像和投影数据[数据集]。癌症成像档案。https://doi.org/10.7937/9npb-2637．

[3]资助国家生物医学成像和生物工程研究所EB017095和EB017185 (Cynthia McCollough, PI)。

[4] Clark, Kenneth, Bruce Vendt, Kirk Smith, John Freymann, Justin Kirby, Paul Koppel, Stephen Moore等人，“癌症成像档案(TCIA):维护和运行一个公共信息库。”数字成像杂志26日,没有。6(2013年12月):1045-57。https://doi.org/10.1007/s10278-013-9622-7．

另请参阅