Семантическая сегментация мультиспектральных изображений с использованием глубокого обучения

Этот пример использует:

Этот пример показывает, как выполнить семантическую сегментацию мультиспектрального изображения с семью каналами, используя U-Net.

Семантическая сегментация включает маркировку каждого пикселя в изображении классом. Одним из вариантов семантической сегментации является отслеживание обезлесения, которое представляет собой изменение лесного покрова с течением времени. Природоохранные учреждения отслеживают обезлесение для оценки и количественной оценки экологического и экологического здоровья области.

Семантическая сегментация на основе глубокого обучения может привести к точному измерению растительного покрова с помощью воздушных фотографий с высоким разрешением. Одной из задач является дифференцирование классов с аналогичными визуальными характеристиками, такими как попытка классифицировать зеленый пиксель как траву, кустарник или дерево. Чтобы повысить точность классификации, некоторые наборы данных содержат мультиспектральные изображения, которые предоставляют дополнительную информацию о каждом пикселе. Например, набор данных Hamlin Beach State Park дополняет цветные изображения тремя ближними инфракрасными каналами, которые обеспечивают более четкое разделение классов.

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

Загрузка данных

Этот пример использует набор мультиспектральных данных высокого разрешения для обучения сети [1]. Набор изображений был запечатлен с помощью беспилотника над государственным парком Хэмлин-Бич, штат Нью-Йорк. Данные содержат маркированные наборы для обучения, валидации и тестирования с 18 метками класса объекта. Размер файла данных составляет ~ 3,0 ГБ.

Загрузите версию MAT-файла набора данных с помощью downloadHamlinBeachMSIData вспомогательная функция. Эта функция присоединена к примеру как вспомогательный файл.

imageDir = tempdir;
url = 'http://www.cis.rit.edu/~rmk6217/rit18_data.mat';
downloadHamlinBeachMSIData(url,imageDir);

Смотрите обучающие данные

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

load(fullfile(imageDir,'rit18_data','rit18_data.mat'));

Исследуйте структуру данных.

whos train_data val_data test_data

  Name            Size                         Bytes  Class     Attributes

  test_data       7x12446x7654            1333663576  uint16              
  train_data      7x9393x5642              741934284  uint16              
  val_data        7x8833x6918              855493716  uint16

Мультиспектральные данные изображения расположены как массивы numChannels-by-width-by-height. Однако в MATLAB ® многоканальные изображения расположены как массивы width-by-height-by-numChannels. Чтобы изменить форму данных так, чтобы каналы находились в третьей размерности, используйте функцию helper, switchChannelsToThirdPlane. Эта функция присоединена к примеру как вспомогательный файл.

train_data = switchChannelsToThirdPlane(train_data);
val_data   = switchChannelsToThirdPlane(val_data);
test_data  = switchChannelsToThirdPlane(test_data);

Подтвердите, что данные имеют правильную структуру.

whos train_data val_data test_data

  Name                Size                     Bytes  Class     Attributes

  test_data       12446x7654x7            1333663576  uint16              
  train_data       9393x5642x7             741934284  uint16              
  val_data         8833x6918x7             855493716  uint16

Цветовыми каналами RGB являются 3-й, 2-й и 1-й каналы изображений. Отобразите цветовой компонент обучающих, валидационных и тестовых изображений в качестве монтажа. Чтобы изображения появлялись ярче на экране, выровняйте их гистограммы с помощью histeq (Image Processing Toolbox) функция.

figure
montage(...
    {histeq(train_data(:,:,[3 2 1])), ...
    histeq(val_data(:,:,[3 2 1])), ...
    histeq(test_data(:,:,[3 2 1]))}, ...
    'BorderSize',10,'BackgroundColor','white')
title('RGB Component of Training Image (Left), Validation Image (Center), and Test Image (Right)')

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

figure
montage(...
    {histeq(train_data(:,:,4)), ...
    histeq(train_data(:,:,5)), ...
    histeq(train_data(:,:,6))}, ...
    'BorderSize',10,'BackgroundColor','white')
title('IR Channels 1 (Left), 2, (Center), and 3 (Right) of Training Image')

Канал 7 является маской, которая указывает на допустимую область сегментации. Отображение маски для изображений для обучения, валидации и тестирования.

figure
montage(...
    {train_data(:,:,7), ...
    val_data(:,:,7), ...
    test_data(:,:,7)}, ...
    'BorderSize',10,'BackgroundColor','white')
title('Mask of Training Image (Left), Validation Image (Center), and Test Image (Right)')

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

disp(classes)

0. Other Class/Image Border      
1. Road Markings                 
2. Tree                          
3. Building                      
4. Vehicle (Car, Truck, or Bus)  
5. Person                        
6. Lifeguard Chair               
7. Picnic Table                  
8. Black Wood Panel              
9. White Wood Panel              
10. Orange Landing Pad           
11. Water Buoy                   
12. Rocks                        
13. Other Vegetation             
14. Grass                        
15. Sand                         
16. Water (Lake)                 
17. Water (Pond)                 
18. Asphalt (Parking Lot/Walkway)

Создайте вектор имен классов.

classNames = [ "RoadMarkings","Tree","Building","Vehicle","Person", ...
               "LifeguardChair","PicnicTable","BlackWoodPanel",...
               "WhiteWoodPanel","OrangeLandingPad","Buoy","Rocks",...
               "LowLevelVegetation","Grass_Lawn","Sand_Beach",...
               "Water_Lake","Water_Pond","Asphalt"];

Наложите метки на обучающее изображение RGB с выравниванием гистограммы. Добавьте шкалу палитры к изображению.

cmap = jet(numel(classNames));
B = labeloverlay(histeq(train_data(:,:,4:6)),train_labels,'Transparency',0.8,'Colormap',cmap);

figure
title('Training Labels')
imshow(B)
N = numel(classNames);
ticks = 1/(N*2):1/N:1;
colorbar('TickLabels',cellstr(classNames),'Ticks',ticks,'TickLength',0,'TickLabelInterpreter','none');
colormap(cmap)

Сохраните обучающие данные как файл MAT, а обучающие метки как файл PNG.

save('train_data.mat','train_data');
imwrite(train_labels,'train_labels.png');

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

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

Начните с хранения обучающих изображений из 'train_data.mat' в imageDatastore. Поскольку формат файла MAT является нестандартным форматом изображения, для чтения данных необходимо использовать средство чтения файлов MAT. Можно использовать вспомогательную программу чтения файлов MAT, matReader, который извлекает первые шесть каналов из обучающих данных и опускает последний канал, содержащий маску. Эта функция присоединена к примеру как вспомогательный файл.

imds = imageDatastore('train_data.mat','FileExtensions','.mat','ReadFcn',@matReader);

Создайте pixelLabelDatastore (Computer Vision Toolbox) для хранения закрашенных фигур меток, содержащих 18 маркированных областей.

pixelLabelIds = 1:18;
pxds = pixelLabelDatastore('train_labels.png',classNames,pixelLabelIds);

Создайте randomPatchExtractionDatastore (Image Processing Toolbox) из datastore изображений и хранилища данных меток пикселей. Каждый мини-пакет содержит 16 закрашенные фигуры размером 256 на 256 пикселей. При каждой итерации эпохи извлекается по тысяче мини-пакетов.

dsTrain = randomPatchExtractionDatastore(imds,pxds,[256,256],'PatchesPerImage',16000);

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

inputBatch = preview(dsTrain);
disp(inputBatch)

        InputImage        ResponsePixelLabelImage
    __________________    _______________________

    {256×256×6 uint16}     {256×256 categorical} 
    {256×256×6 uint16}     {256×256 categorical} 
    {256×256×6 uint16}     {256×256 categorical} 
    {256×256×6 uint16}     {256×256 categorical} 
    {256×256×6 uint16}     {256×256 categorical} 
    {256×256×6 uint16}     {256×256 categorical} 
    {256×256×6 uint16}     {256×256 categorical} 
    {256×256×6 uint16}     {256×256 categorical}

Создание слоев сети U-Net

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

Этот пример изменяет U-Net, чтобы использовать заполнение нулями в свертках, так что вход и выход в свертки имеют одинаковый размер. Используйте функцию helper, createUnet, чтобы создать U-Net с несколькими предварительно выбранными гиперпараметрами. Эта функция присоединена к примеру как вспомогательный файл.

inputTileSize = [256,256,6];
lgraph = createUnet(inputTileSize);
disp(lgraph.Layers)

  58x1 Layer array with layers:

     1   'ImageInputLayer'                        Image Input                  256x256x6 images with 'zerocenter' normalization
     2   'Encoder-Section-1-Conv-1'               Convolution                  64 3x3x6 convolutions with stride [1  1] and padding [1  1  1  1]
     3   'Encoder-Section-1-ReLU-1'               ReLU                         ReLU
     4   'Encoder-Section-1-Conv-2'               Convolution                  64 3x3x64 convolutions with stride [1  1] and padding [1  1  1  1]
     5   'Encoder-Section-1-ReLU-2'               ReLU                         ReLU
     6   'Encoder-Section-1-MaxPool'              Max Pooling                  2x2 max pooling with stride [2  2] and padding [0  0  0  0]
     7   'Encoder-Section-2-Conv-1'               Convolution                  128 3x3x64 convolutions with stride [1  1] and padding [1  1  1  1]
     8   'Encoder-Section-2-ReLU-1'               ReLU                         ReLU
     9   'Encoder-Section-2-Conv-2'               Convolution                  128 3x3x128 convolutions with stride [1  1] and padding [1  1  1  1]
    10   'Encoder-Section-2-ReLU-2'               ReLU                         ReLU
    11   'Encoder-Section-2-MaxPool'              Max Pooling                  2x2 max pooling with stride [2  2] and padding [0  0  0  0]
    12   'Encoder-Section-3-Conv-1'               Convolution                  256 3x3x128 convolutions with stride [1  1] and padding [1  1  1  1]
    13   'Encoder-Section-3-ReLU-1'               ReLU                         ReLU
    14   'Encoder-Section-3-Conv-2'               Convolution                  256 3x3x256 convolutions with stride [1  1] and padding [1  1  1  1]
    15   'Encoder-Section-3-ReLU-2'               ReLU                         ReLU
    16   'Encoder-Section-3-MaxPool'              Max Pooling                  2x2 max pooling with stride [2  2] and padding [0  0  0  0]
    17   'Encoder-Section-4-Conv-1'               Convolution                  512 3x3x256 convolutions with stride [1  1] and padding [1  1  1  1]
    18   'Encoder-Section-4-ReLU-1'               ReLU                         ReLU
    19   'Encoder-Section-4-Conv-2'               Convolution                  512 3x3x512 convolutions with stride [1  1] and padding [1  1  1  1]
    20   'Encoder-Section-4-ReLU-2'               ReLU                         ReLU
    21   'Encoder-Section-4-DropOut'              Dropout                      50% dropout
    22   'Encoder-Section-4-MaxPool'              Max Pooling                  2x2 max pooling with stride [2  2] and padding [0  0  0  0]
    23   'Mid-Conv-1'                             Convolution                  1024 3x3x512 convolutions with stride [1  1] and padding [1  1  1  1]
    24   'Mid-ReLU-1'                             ReLU                         ReLU
    25   'Mid-Conv-2'                             Convolution                  1024 3x3x1024 convolutions with stride [1  1] and padding [1  1  1  1]
    26   'Mid-ReLU-2'                             ReLU                         ReLU
    27   'Mid-DropOut'                            Dropout                      50% dropout
    28   'Decoder-Section-1-UpConv'               Transposed Convolution       512 2x2x1024 transposed convolutions with stride [2  2] and cropping [0  0  0  0]
    29   'Decoder-Section-1-UpReLU'               ReLU                         ReLU
    30   'Decoder-Section-1-DepthConcatenation'   Depth concatenation          Depth concatenation of 2 inputs
    31   'Decoder-Section-1-Conv-1'               Convolution                  512 3x3x1024 convolutions with stride [1  1] and padding [1  1  1  1]
    32   'Decoder-Section-1-ReLU-1'               ReLU                         ReLU
    33   'Decoder-Section-1-Conv-2'               Convolution                  512 3x3x512 convolutions with stride [1  1] and padding [1  1  1  1]
    34   'Decoder-Section-1-ReLU-2'               ReLU                         ReLU
    35   'Decoder-Section-2-UpConv'               Transposed Convolution       256 2x2x512 transposed convolutions with stride [2  2] and cropping [0  0  0  0]
    36   'Decoder-Section-2-UpReLU'               ReLU                         ReLU
    37   'Decoder-Section-2-DepthConcatenation'   Depth concatenation          Depth concatenation of 2 inputs
    38   'Decoder-Section-2-Conv-1'               Convolution                  256 3x3x512 convolutions with stride [1  1] and padding [1  1  1  1]
    39   'Decoder-Section-2-ReLU-1'               ReLU                         ReLU
    40   'Decoder-Section-2-Conv-2'               Convolution                  256 3x3x256 convolutions with stride [1  1] and padding [1  1  1  1]
    41   'Decoder-Section-2-ReLU-2'               ReLU                         ReLU
    42   'Decoder-Section-3-UpConv'               Transposed Convolution       128 2x2x256 transposed convolutions with stride [2  2] and cropping [0  0  0  0]
    43   'Decoder-Section-3-UpReLU'               ReLU                         ReLU
    44   'Decoder-Section-3-DepthConcatenation'   Depth concatenation          Depth concatenation of 2 inputs
    45   'Decoder-Section-3-Conv-1'               Convolution                  128 3x3x256 convolutions with stride [1  1] and padding [1  1  1  1]
    46   'Decoder-Section-3-ReLU-1'               ReLU                         ReLU
    47   'Decoder-Section-3-Conv-2'               Convolution                  128 3x3x128 convolutions with stride [1  1] and padding [1  1  1  1]
    48   'Decoder-Section-3-ReLU-2'               ReLU                         ReLU
    49   'Decoder-Section-4-UpConv'               Transposed Convolution       64 2x2x128 transposed convolutions with stride [2  2] and cropping [0  0  0  0]
    50   'Decoder-Section-4-UpReLU'               ReLU                         ReLU
    51   'Decoder-Section-4-DepthConcatenation'   Depth concatenation          Depth concatenation of 2 inputs
    52   'Decoder-Section-4-Conv-1'               Convolution                  64 3x3x128 convolutions with stride [1  1] and padding [1  1  1  1]
    53   'Decoder-Section-4-ReLU-1'               ReLU                         ReLU
    54   'Decoder-Section-4-Conv-2'               Convolution                  64 3x3x64 convolutions with stride [1  1] and padding [1  1  1  1]
    55   'Decoder-Section-4-ReLU-2'               ReLU                         ReLU
    56   'Final-ConvolutionLayer'                 Convolution                  18 1x1x64 convolutions with stride [1  1] and padding [0  0  0  0]
    57   'Softmax-Layer'                          Softmax                      softmax
    58   'Segmentation-Layer'                     Pixel Classification Layer   Cross-entropy loss

Выбор опций обучения

Обучите сеть с помощью стохастического градиентного спуска с оптимизацией импульса (SGDM). Задайте установки гиперпараметров для SGDM при помощи trainingOptions функция.

Обучение глубокой сети занимает много времени. Ускорите обучение, задав высокую скорость обучения. Однако это может заставить градиенты сети взрываться или расти бесконтрольно, препятствуя успешному обучению сети. Чтобы сохранить градиенты в значимой области значений, включите усечение градиента путем определения 'GradientThreshold' как 0.05, и задайте 'GradientThresholdMethod' для использования L2-norm градиентов.

initialLearningRate = 0.05;
maxEpochs = 150;
minibatchSize = 16;
l2reg = 0.0001;

options = trainingOptions('sgdm',...
    'InitialLearnRate',initialLearningRate, ...
    'Momentum',0.9,...
    'L2Regularization',l2reg,...
    'MaxEpochs',maxEpochs,...
    'MiniBatchSize',minibatchSize,...
    'LearnRateSchedule','piecewise',...    
    'Shuffle','every-epoch',...
    'GradientThresholdMethod','l2norm',...
    'GradientThreshold',0.05, ...
    'Plots','training-progress', ...
    'VerboseFrequency',20);

Обучите сеть

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

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

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

doTraining = false; 
if doTraining
    [net,info] = trainNetwork(dsTrain,lgraph,options);
    modelDateTime = string(datetime('now','Format',"yyyy-MM-dd-HH-mm-ss"));
    save(strcat("multispectralUnet-",modelDateTime,"-Epoch-",num2str(maxEpochs),".mat"),'net');

else 
    trainedUnet_url = 'https://www.mathworks.com/supportfiles/vision/data/multispectralUnet.mat';
    downloadTrainedUnet(trainedUnet_url,imageDir);
    load(fullfile(imageDir,'trainedUnet','multispectralUnet.mat'));
end

Теперь можно использовать U-Net, чтобы семантически сегментировать мультиспектральное изображение.

Предсказание результатов тестовых данных

Чтобы выполнить прямой проход по обученной сети, используйте функцию helper, segmentImage, с набором данных валидации. Эта функция присоединена к примеру как вспомогательный файл. segmentImage выполняет сегментацию на закрашенных фигурах к изображениям с помощью semanticseg (Computer Vision Toolbox) функция.

predictPatchSize = [1024 1024];
segmentedImage = segmentImage(val_data,net,predictPatchSize);

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

segmentedImage = uint8(val_data(:,:,7)~=0) .* segmentedImage;

figure
imshow(segmentedImage,[])
title('Segmented Image')

Выход семантической сегментации зашумлен. Выполните постобработку изображений, чтобы удалить шум и шальные пиксели. Используйте medfilt2 (Image Processing Toolbox) функция для удаления шума от соли и перца из сегментации. Визуализируйте сегментированное изображение с удаленным шумом.

segmentedImage = medfilt2(segmentedImage,[7,7]);
imshow(segmentedImage,[]);
title('Segmented Image  with Noise Removed')

Наложите сегментированное изображение на изображение валидации RGB с выравниванием гистограммы.

B = labeloverlay(histeq(val_data(:,:,[3 2 1])),segmentedImage,'Transparency',0.8,'Colormap',cmap);

figure
imshow(B)
title('Labeled Validation Image')
colorbar('TickLabels',cellstr(classNames),'Ticks',ticks,'TickLength',0,'TickLabelInterpreter','none');
colormap(cmap)

Сохраните сегментированное изображение и основную истину метки как файлы PNG. Они будут использоваться для вычисления метрик точности.

imwrite(segmentedImage,'results.png');
imwrite(val_labels,'gtruth.png');

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

Создайте pixelLabelDatastore (Computer Vision Toolbox) для результатов сегментации и наземных меток истинности.

pxdsResults = pixelLabelDatastore('results.png',classNames,pixelLabelIds);
pxdsTruth = pixelLabelDatastore('gtruth.png',classNames,pixelLabelIds);

Измерьте глобальную точность семантической сегментации с помощью evaluateSemanticSegmentation (Computer Vision Toolbox) функция.

ssm = evaluateSemanticSegmentation(pxdsResults,pxdsTruth,'Metrics','global-accuracy');

Evaluating semantic segmentation results
----------------------------------------
* Selected metrics: global accuracy.
* Processed 1 images.
* Finalizing... Done.
* Data set metrics:

    GlobalAccuracy
    ______________

       0.90698

Глобальный счет точности указывает, что чуть более 90% пикселей классифицированы правильно.

Вычисление объема растительного покрова

Конечной целью этого примера является вычисление степени растительного покрова на мультиспектральном изображении.

Найдите количество пикселей, маркированных растительностью. Идентификаторы меток 2 («Деревья»), 13 («LowLevelVegetation») и 14 («Grass _ Lawn») являются классами растительности. Также найдите общее количество допустимых пикселей путем суммирования пикселей в информация только для чтения маскировочного изображения.

vegetationClassIds = uint8([2,13,14]);
vegetationPixels = ismember(segmentedImage(:),vegetationClassIds);
validPixels = (segmentedImage~=0);

numVegetationPixels = sum(vegetationPixels(:));
numValidPixels = sum(validPixels(:));

Вычислите процент растительного покрова путем деления количества пикселей растительности на количество допустимых пикселей.

percentVegetationCover = (numVegetationPixels/numValidPixels)*100;
fprintf('The percentage of vegetation cover is %3.2f%%.',percentVegetationCover);

The percentage of vegetation cover is 51.72%.

Ссылки

[1] Kemker, R., C. Salvaggio, and C. Kanan. «Мультиспектральный набор данных высокого разрешения для семантической сегментации». CoRR, abs/1703.01918. 2017.

[2] Роннебергер, О., П. Фишер и Т. Брокс. U-Net: Сверточные сети для сегментации биомедицинских изображений. CoRR, abs/1505.04597. 2015.

См. также

Внешние веб-сайты

https://github.com/rmkemker/RIT-18

Документация