MATLAB环境下基于深度学习的多光谱图像语义分割

深度学习，一个大号/深层的，现代的，黑箱的，信号/图像处理器，本文程序运行环境为MATLAB R2018A。

本文简要讲解如何训练U-Net卷积神经网络对7个通道的多光谱图像进行语义分割，7个通道包括3个颜色通道、3个近红外通道和一个掩模通道。需要使用具有计算能力3.0或更高版本的支持CUDA的NVIDIA?GPU（需要并行计算工具箱）。

代码链接

https://mianbaoduo.com/o/bread/mbd-Y5mYlJ5x

引言

语义分割使用一个类别来标记图像中的每个像素，一个很典型的应用是跟踪森林覆盖度随时间的变化，以评估和量化一个地区的环境和生态健康状况。

基于深度学习的语义分割可以从高分辨率的航空照片中精确地测量植被覆盖度。其中的一个难点是区分具有相似视觉特征的类别，例如将一个绿色像素分类为草、灌木或树。为了提高分类精度，数据集包含多光谱图像，提供关于每个像素的额外信息。例如，哈姆林海滩州立公园的数据集用近红外通道补充了彩色图像，提供了更清晰的类别划分。

本文简单讲解如何使用基于深度学习的语义分割技术从一组多光谱图像中计算一个区域的植被覆盖百分比。

本文使用高分辨率的多光谱数据集来训练深层网络。数据集是用无人机在纽约州哈姆林海滩州立公园上空拍摄的，包含标记的训练、验证和测试集，带有18个对象类别标签，数据为MAT文件，大小约为3.0 GB。

检查训练数据

将数据集加载到MATLAB工作区

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 的数组，因此需要reshape数据，使通道处于第三个维度中。

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颜色通道是第4、5、6个图像通道，展示训练、验证和测试图像的颜色分量，为了使图像在屏幕上看起来更亮，使用histeq平衡直方图

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

显示训练数据的前3个直方图均衡通道，这些通道对应于近红外波段，并根据它们的热特征突出显示图像的不同组成部分。例如，靠近第二通道图像中心的树比其他两个通道的树显示出更多的细节。

figure
montage(...
 {histeq(train_data(:,:,1)), ...
 histeq(train_data(:,:,2)), ...
 histeq(train_data(:,:,3))}, ...
 '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个类别中的一个，同时获取对应id的类别的列表。

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');

定义Mini-Batch Datastore

使用一个mini-batch datastore将训练数据输入网络，本例定义mini-batch datastore的自定义实现，即PixelLabelImagePatchDatastore，作为一种方便生成增强图像的方法，用于训练语义分割网络。

将“train_data.mat”的训练图像存储在 imageDatastore中，由于MAT文件格式是一种非标准的图像格式，因此必须使用MAT文件阅读器读取图像数据，本例使用MAT文件阅读器matReader，它从训练数据中提取前六个通道，并省略包含掩码通道的最后一个通道。

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

创建一个 pixelLabelDatastore 存储包含18个标记区域的“标签块”

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

PixelLabelImagePatchDatastore函数从训练图像和标签中提取块(patches)对，此外还对训练数据进行数据增强，以增加训练数据集的大小，并使训练更具鲁棒性。为了模拟无人机收集的空中景观数据，PixelLabelImagePatchDatastore 通过添加随机旋转、反射和缩放来增加块(patches)对。每个mini-batch包含16个大小为256*256像素的块（patches），所有的块（patches）都是从图像中的随机位置中提取出来的。

ds = PixelLabelImagePatchDatastore(imds,pxds,'PatchSize',256,...
 'MiniBatchSize',16,...
 'BatchesPerImage',1000) 

设置U-Net网络

本例使用U-Net网络的一个变体，在U-Net网络中，初始一系列的卷积层与最大池化层穿插起来，依次降低了输入图像的分辨率。这些层之后是一系列的卷积层，中间穿插着上采样算子，依次提高了输入图像的分辨率。本例修改了U-Net，在卷积中使用零填充，以便使卷积的输入和输出具有相同的大小。使用函数createUnet创建一个具有预选超参数的U-Net。

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

训练参数设置

使用动量随机梯度下降（SGDM）优化算法来训练网络，通过使用 trainingOptions 函数指定SDGM的超参数设置。

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

训练网络

在配置完训练参数和 mini-batch datastore后，使用trainNetwork函数对U-Net网络进行训练。

注意：在NVIDIA（TM）Titan X上进行训练大约需要20个小时

 [net,info] = trainNetwork(datasource,lgraph,options); 

现在可以使用U-Net对多光谱图像进行语义分割。

模型测试

使用辅助函数segmentImage分割图像

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

为了只提取分割的有效部分，将分割后的图像乘以验证数据的掩码通道。

segmentedImage = uint8(val_data(:,:,7)~=0) .* segmentedImage;
figure
imshow(segmentedImage,[])
title('Segmented Image')

语义分割的输出是有噪声的，因此要执行图像后处理，以去除噪声和杂散像素。使用medfilt2函数去除椒盐噪声，并可视化分割后的图像。

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

将分割后的图像叠加在直方图均衡的RGB验证图像上

B = labeloverlay(histeq(val_data(:,:,4:6)),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 数据存储

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

利用 evaluateSemanticSegmentation 函数计算语义分割的全局精度。

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

Evaluating semantic segmentation results
----------------------------------------
* Selected metrics: global accuracy.
* Processing 1 images...
[==================================================] 100%
Elapsed time: 00:00:25
Estimated time remaining: 00:00:00
* Finalizing... Done.
* Data set metrics:

GlobalAccuracy
______________

0.90698

全局精度得分表明，超过90%的像素被正确分类。

计算植被覆盖度的范围

本例的最终目标是计算多光谱图像中的植被覆盖程度。

找到标记为植被的像素数。标签2（“树木”）、13（“低水平植被”）和14（“草坪”）是植被类别。同时，通过将掩模图像的ROI中的像素相加，也可以找到有效像素的总数。

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. "High-Resolution Multispectral Dataset for Semantic Segmentation." CoRR, abs/1703.01918. 2017.

[2] Ronneberger, O., P. Fischer, and T. Brox. "U-Net: Convolutional Networks for Biomedical Image Segmentation." CoRR, abs/1505.04597. 2015.