Загрузите набор данных LDCT

Этот пример использует данные из Низкого CT Дозы Главная проблема [2, 3, 4]. Данные включают пары изображений CT регулярной дозы, и симулированные изображения CT низкой дозы для 99 главных сканов (пометил N для neuro), 100 сканов грудной клетки (пометил C для груди), и 100 сканов живота (пометил L для печени).

Установите dataDir как желаемое местоположение набора данных. Данные для этого примера требуют 52 Гбайт памяти.

dataDir = fullfile(tempdir,"LDCT","LDCT-and-Projection-data");

Чтобы загрузить данные, перейдите к веб-сайту Архива Обработки изображений Рака. Этот пример использует ony два терпеливых скана из груди. Загрузите файлы "C081" и "C120" грудной клетки с "Изображений (DICOM, 952 Гбайт)" набор данных с помощью Ретривера Данных 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 позволять считать данные.

exts = {'.dcm'};
readFcn = @(x)rescale(dicomread(x));
imdsLDTrain = imageDatastore(fullfile(scanDirTrain,"1.000000-Low Dose Images-71581"), ...
    FileExtensions=exts,ReadFcn=readFcn);
imdsHDTrain = imageDatastore(fullfile(scanDirTrain,"1.000000-Full dose images-34601"), ...
    FileExtensions=exts,ReadFcn=readFcn);
imdsLDTest = imageDatastore(fullfile(scanDirTest,"1.000000-Low Dose Images-32837"), ...
    FileExtensions=exts,ReadFcn=readFcn);
imdsHDTest = imageDatastore(fullfile(scanDirTest,"1.000000-Full dose images-95670"), ...
    FileExtensions=exts,ReadFcn=readFcn);

Предварительно просмотрите учебное изображение от низкой дозы и обучающих наборов данных CT большей дозы.

lowDose = preview(imdsLDTrain);
highDose = preview(imdsHDTrain);
montage({lowDose,highDose})

Предварительно обработайте и увеличьте обучающие данные

Задайте входной размер изображений для входных и выходных изображений.

inputSize = [256,256,1];

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

augmentDataForLD2HDCT функция выполняет эти операции:

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

Набор данных LDCT обеспечивает пары низкой дозы и изображений CT большей дозы. Однако МОДУЛЬНАЯ архитектура требует непарных данных для безнадзорного изучения. Этот пример симулирует непарные данные об обучении и валидации путем перестановки данных в каждой итерации учебного цикла.

Обработайте данные об обучении и валидации в пакетном режиме во время обучения

Этот пример использует пользовательский учебный цикл. minibatchqueue объект полезен для управления мини-пакетная обработка наблюдений в пользовательских учебных циклах. minibatchqueue возразите также бросает данные к dlarray объект, который включает автоматическое дифференцирование в применении глубокого обучения.

Задайте мини-пакетный формат экстракции данных как SSCB (пространственный, пространственный, канал, пакет). Установите DispatchInBackground аргумент значения имени как boolean, возвращенный canUseGPU. Если поддерживаемый графический процессор доступен для расчета, то minibatchqueue объект предварительно обрабатывает мини-пакеты в фоновом режиме в параллельном пуле во время обучения.

miniBatchSize = 1;
mbqLDTrain = minibatchqueue(imdsLDTrain,MiniBatchSize=miniBatchSize, ...
    MiniBatchFormat="SSCB",DispatchInBackground=canUseGPU);
mbqHDTrain = minibatchqueue(imdsHDTrain,MiniBatchSize=miniBatchSize, ...
    MiniBatchFormat="SSCB",DispatchInBackground=canUseGPU);

Создайте сеть генератора

МОДУЛЬ состоит из одного генератора и двух различителей. Генератор выполняет перевод от изображения к изображению от низкой дозы до большей дозы. Различители являются сетями PatchGAN, которые возвращают мудрую закрашенной фигурой вероятность, что входные данные действительны или сгенерированы. Один различитель различает действительные и сгенерированные изображения низкой дозы, и другой различитель различает действительные и сгенерированные изображения большей дозы.

Создайте МОДУЛЬНУЮ сеть генератора использование unitGenerator (Image Processing Toolbox) функция. Входные и выходные разделы энкодера генератора каждый состоит из двух блоков субдискретизации и пяти остаточных блоков. Разделы энкодера совместно используют два из пяти остаточных блоков. Аналогично, входные и выходные разделы декодера генератора, каждый состоит из двух блоков субдискретизации и пяти остаточных блоков и разделов декодера, совместно используют два из пяти остаточных блоков.

gen = unitGenerator(inputSize);

Визуализируйте сеть генератора.

analyzeNetwork(gen)

Создайте сети различителя

Существует две сети различителя, один для каждой из областей изображения (CT низкой дозы и CT большей дозы). Создайте различители для входных и выходных областей с помощью patchGANDiscriminator (Image Processing Toolbox) функция.

discLD = patchGANDiscriminator(inputSize,NumDownsamplingBlocks=4,FilterSize=3, ...
    ConvolutionWeightsInitializer="narrow-normal",NormalizationLayer="none");
discHD = patchGANDiscriminator(inputSize,"NumDownsamplingBlocks",4,FilterSize=3, ...
    ConvolutionWeightsInitializer="narrow-normal",NormalizationLayer="none");

Визуализируйте сети различителя.

analyzeNetwork(discLD);
analyzeNetwork(discHD);

Градиенты модели Define и функции потерь

modelGradientDisc и modelGradientGen функции помощника вычисляют градиенты и потери для различителей и генератора, соответственно. Эти функции заданы в разделе Supporting Functions этого примера.

Цель каждого различителя состоит в том, чтобы правильно различать действительные изображения (1) и переведенные изображения (0) для изображений в его области. Каждый различитель имеет одну функцию потерь.

Цель генератора состоит в том, чтобы сгенерировать переведенные изображения, которые различители классифицируют как действительные. Потеря генератора является взвешенной суммой пяти типов потерь: потеря самореконструкции, потеря непротиворечивости цикла, скрытая потеря KL, цикл скрытая потеря KL и соперничающая потеря.

Задайте весовые коэффициенты за различные потери.

lossWeights.selfReconLossWeight = 10;
lossWeights.hiddenKLLossWeight = 0.01;
lossWeights.cycleConsisLossWeight = 10;
lossWeights.cycleHiddenKLLossWeight = 0.01;
lossWeights.advLossWeight = 1;
lossWeights.discLossWeight = 0.5;

Задайте опции обучения

Задайте опции для оптимизации Адама. Обучите сеть в течение 26 эпох.

numEpochs = 26;

Задайте идентичные опции для сетей различителя и генератора.

learnRate = 0.0001;
gradDecay = 0.5;
sqGradDecay = 0.999;
weightDecay = 0.0001;

genAvgGradient = [];
genAvgGradientSq = [];
discLDAvgGradient = [];
discLDAvgGradientSq = [];
discHDAvgGradient = [];
discHDAvgGradientSq = [];

Обучите модель или загрузите предварительно обученную МОДУЛЬНУЮ сеть

По умолчанию пример загружает предварительно обученную версию МОДУЛЬНОГО генератора для набора данных CT NIH-AAPM-Mayo Low-Dose при помощи функции помощника downloadTrainedLD2HDCTUNITNet. Функция помощника присоединена к примеру как к вспомогательному файлу. Предварительно обученная сеть позволяет вам запустить целый пример, не ожидая обучения завершиться.

Чтобы обучить сеть, установите doTraining переменная в следующем коде к true. Обучите модель в пользовательском учебном цикле. Для каждой итерации:

Обучайтесь на графическом процессоре, если вы доступны. Используя графический процессор требует Parallel Computing Toolbox™, и CUDA® включил NVIDIA® графический процессор. Для получения дополнительной информации смотрите Поддержку графического процессора Релизом (Parallel Computing Toolbox). Обучение занимает приблизительно 58 часов на Титане NVIDIA RTX.

doTraining = false;
if doTraining

    % Create a figure to show the results
    figure("Units","Normalized");
    for iPlot = 1:4
        ax(iPlot) = subplot(2,2,iPlot);
    end
    
    iteration = 0;

    % Loop over epochs
    for epoch = 1:numEpochs
        
        % Shuffle data every epoch
        reset(mbqLDTrain);
        shuffle(mbqLDTrain);
        reset(mbqHDTrain);
        shuffle(mbqHDTrain);
        
        % Run the loop until all the images in the mini-batch queue
        % mbqLDTrain are processed
        while hasdata(mbqLDTrain)
            iteration = iteration + 1;
            
            % Read data from the low-dose domain
            imLowDose = next(mbqLDTrain); 
             
            % Read data from the high-dose domain
            if hasdata(mbqHDTrain) == 0
                reset(mbqHDTrain);
                shuffle(mbqHDTrain);
            end
            imHighDose = next(mbqHDTrain);
    
            % Calculate discriminator gradients and losses
            [discLDGrads,discHDGrads,discLDLoss,discHDLoss] = dlfeval(@modelGradientDisc, ...
                gen,discLD,discHD,imLowDose,imHighDose,lossWeights.discLossWeight);
            
            % Apply weight decay regularization on low-dose discriminator gradients
            discLDGrads = dlupdate(@(g,w) g+weightDecay*w,discLDGrads,discLD.Learnables);
            
            % Update parameters of low-dose discriminator
            [discLD,discLDAvgGradient,discLDAvgGradientSq] = adamupdate(discLD,discLDGrads, ...
                discLDAvgGradient,discLDAvgGradientSq,iteration,learnRate,gradDecay,sqGradDecay);  
            
            % Apply weight decay regularization on high-dose discriminator gradients
            discHDGrads = dlupdate(@(g,w) g+weightDecay*w,discHDGrads,discHD.Learnables);
            
            % Update parameters of high-dose discriminator
            [discHD,discHDAvgGradient,discHDAvgGradientSq] = adamupdate(discHD,discHDGrads, ...
                discHDAvgGradient,discHDAvgGradientSq,iteration,learnRate,gradDecay,sqGradDecay);
            
            % Calculate generator gradient and loss
            [genGrad,genLoss,images] = dlfeval(@modelGradientGen,gen,discLD,discHD,imLowDose,imHighDose,lossWeights);
            
            % Apply weight decay regularization on generator gradients
            genGrad = dlupdate(@(g,w) g+weightDecay*w,genGrad,gen.Learnables);
            
            % Update parameters of generator
            [gen,genAvgGradient,genAvgGradientSq] = adamupdate(gen,genGrad,genAvgGradient, ...
                genAvgGradientSq,iteration,learnRate,gradDecay,sqGradDecay);
        end
        
        % Display the results
        updateTrainingPlotLowDoseToHighDose(ax,images{:});
    end
    
    % Save the trained network
    modelDateTime = string(datetime("now",Format="yyyy-MM-dd-HH-mm-ss"));
    save(strcat("trainedLowDoseHighDoseUNITGeneratorNet-",modelDateTime,"-Epoch-",num2str(numEpochs),".mat"),'gen');
    
else
    net_url = "https://ssd.mathworks.com/supportfiles/vision/data/trainedLowDoseHighDoseUNITGeneratorNet.zip";
    downloadTrainedLD2HDCTUNITNet(net_url,dataDir);
    load(fullfile(dataDir,"trainedLowDoseHighDoseUNITGeneratorNet.mat"));
end

Сгенерируйте изображение большей дозы Используя, обучил сеть

Считайте и отобразите изображение от datastore тестовых изображений низкой дозы.

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

Преобразуйте изображение в тип данных single. Перемасштабируйте данные изображения к области значений [-1, 1] как ожидалось последним слоем сети генератора.

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

Создайте dlarray возразите что входные данные против генератора. Если поддерживаемый графический процессор доступен для расчета, то выполните вывод на графическом процессоре путем преобразования данных в gpuArray объект.

dlLowDoseImage = dlarray(imLowDoseTestRescaled,'SSCB');    
if canUseGPU
    dlLowDoseImage = gpuArray(dlLowDoseImage);
end

Переведите входное изображение низкой дозы в область большей дозы использование unitPredict (Image Processing Toolbox) функция. Сгенерированное изображение имеет пиксельные значения в области значений [-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(['Low-dose Test Image ',num2str(idxToTest),' with Generated High-dose Image and Ground Truth High-dose Image'])

Вспомогательные Функции

Функция градиентов модели

modelGradientGen функция помощника вычисляет градиенты и потерю для генератора.

function [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");
    
    % Calculate different losses
    selfReconLoss = computeReconLoss(imLD,imLD2LD) + computeReconLoss(imHD,imHD2HD);
    hiddenKLLoss = computeKLLoss(hidden);
    cycleReconLoss = computeReconLoss(imLD,imLD2HD2LD) + computeReconLoss(imHD,imHD2LD2HD);
    cycleHiddenKLLoss = computeKLLoss(cycle_hidden);
    
    outA = forward(discLD,imHD2LD);
    outB = forward(discHD,imLD2HD);
    advLoss = computeAdvLoss(outA) + computeAdvLoss(outB);
    
    % Calculate the total loss of generator as a weighted sum of five losses
    genTotalLoss = ...
        selfReconLoss*lossWeights.selfReconLossWeight + ...
        hiddenKLLoss*lossWeights.hiddenKLLossWeight + ...
        cycleReconLoss*lossWeights.cycleConsisLossWeight + ...
        cycleHiddenKLLoss*lossWeights.cycleHiddenKLLossWeight + ...
        advLoss*lossWeights.advLossWeight;
    
    % Update the parameters of generator
    genGrad = dlgradient(genTotalLoss,gen.Learnables); 
    
    % Convert the data type from dlarray to single
    genLoss = extractdata(genTotalLoss);
    images = {imLD,imLD2HD,imHD,imHD2LD};
end

modelGradientDisc функция помощника вычисляет градиенты и потерю для этих двух различителей.

function [discLDGrads,discHDGrads,discLDLoss,discHDLoss] = modelGradientDisc(gen, ...
    discLD,discHD,imRealLD,imRealHD,discLossWeight)

    [~,imFakeLD,imFakeHD,~] = forward(gen,imRealLD,imRealHD);
    
    % Calculate loss of the discriminator for low-dose images
    outRealLD = forward(discLD,imRealLD); 
    outFakeLD = forward(discLD,imFakeLD);
    discLDLoss = discLossWeight*computeDiscLoss(outRealLD,outFakeLD);
    
    % Update parameters of the discriminator for low-dose images
    discLDGrads = dlgradient(discLDLoss,discLD.Learnables); 
    
    % Calculate loss of the discriminator for high-dose images
    outRealHD = forward(discHD,imRealHD); 
    outFakeHD = forward(discHD,imFakeHD);
    discHDLoss = discLossWeight*computeDiscLoss(outRealHD,outFakeHD);
    
    % Update parameters of the discriminator for high-dose images
    discHDGrads = dlgradient(discHDLoss,discHD.Learnables);
    
    % Convert the data type from dlarray to single
    discLDLoss = extractdata(discLDLoss);
    discHDLoss = extractdata(discHDLoss);
end

Функции потерь

computeDiscLoss функция помощника вычисляет потерю различителя. Каждая потеря различителя является суммой двух компонентов:

function discLoss = computeDiscLoss(Yreal,Ytranslated)
    discLoss = mean(((1-Yreal).^2),"all") + ...
               mean(((0-Ytranslated).^2),"all");
end

computeAdvLoss функция помощника вычисляет соперничающую потерю для генератора. Соперничающая потеря является различием в квадрате между вектором из единиц и предсказаниями различителя на переведенном изображении.

function advLoss = computeAdvLoss(Ytranslated)
    advLoss = mean(((Ytranslated-1).^2),"all");
end

computeReconLoss функция помощника вычисляет потерю самореконструкции и потерю непротиворечивости цикла для генератора. Сам потеря реконструкции ${\mathit{L}}^1$ расстояние между входными изображениями и их самовосстановленными версиями. Потеря непротиворечивости цикла ${\mathit{L}}^1$ расстояние между входными изображениями и их восстановленными циклом версиями.

function reconLoss = computeReconLoss(Yreal,Yrecon)
    reconLoss = mean(abs(Yreal-Yrecon),"all");
end

computeKLLoss функция помощника вычисляет скрытую потерю KL и скрытую от цикла потерю KL для генератора. Скрытая потеря KL является различием в квадрате между нулевым вектором и 'encoderSharedBlock'активация для потока самореконструкции. Скрытая от цикла потеря KL является различием в квадрате между нулевым вектором и 'encoderSharedBlock'активация для потока реконструкции цикла.

function klLoss = computeKLLoss(hidden)
    klLoss = mean(abs(hidden.^2),"all");
end

Ссылки

[1] Лю, Мин-Юй, Томас Бреуель и Ян Коц, "Безнадзорные сети перевода от изображения к изображению". В Усовершенствованиях в Нейронных Системах обработки информации, 2017. https://arxiv.org/pdf/1703.00848.pdf.

[2] МакКалоу, C.H., Чен, B., Холмс, D., III, Дуань, X., Ю, Z., Ю, L., Лэн, S., Флетчер, J. (2020). Данные из Низких Данных об Изображении и Проекции CT Дозы [Набор данных]. Архив Обработки изображений Рака. https://doi.org/10.7937/9npb-2637.

[3] Предоставления EB017095 и EB017185 (Синтия МакКалоу, PI) от национального института биомедицинской обработки изображений и биоинженерии.

[4] Кларк, Кеннет, Брюс Вендт, Кирк Смит, Джон Фреиманн, Джастин Кирби, Пол Коппель, Стивен Мур, и др. “Архив обработки изображений рака (TCIA): Поддержание и Работа Репозиторием Общедоступной информации”. Журнал Цифровой Обработки изображений 26, № 6 (декабрь 2013): 1045–57. https://doi.org/10.1007/s10278-013-9622-7.