Main Content

本页的翻译已过时。点击此处可查看最新英文版本。

使用深度学习对多光谱图像进行语义分割

此示例说明如何使用 U-Net 对包含七个通道的多光谱图像执行语义分割。

语义分割涉及用类标注图像中的每个像素。语义分割的一个应用是跟踪森林采伐,即森林植被随时间的变化。环保机构需要跟踪森林采伐,以评估和量化某个地区的环境和生态健康。

基于深度学习的语义分割可以通过高分辨率航拍照片精确测量植被覆盖度。一个挑战是区分具有相似视觉特性的类,例如尝试将绿色像素分类为草、灌木或树。为了提高分类准确度,一些数据集包含多光谱图像,这些图像提供关于每个像素的附加信息。例如,Hamlin Beach 国家公园数据集用三个近红外通道为彩色图像补充信息,以提供更清晰的分类。

此示例说明如何使用基于深度学习的语义分割方法,根据一组多光谱图像计算某个区域的植被覆盖率。

下载数据

此示例使用高分辨率多光谱数据集来训练网络 [1]。图像集是用无人机在纽约 Hamlin Beach 州立公园拍摄的。数据包含已标注的训练集、验证集和测试集,有 18 个对象类标签。数据文件的大小约为 3.0 GB。

使用 downloadHamlinBeachMSIData 辅助函数下载数据集的 MAT 文件版本。此函数作为支持文件包含在本示例中。

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              

多光谱图像数据排列为通道数×宽度×高度数组。但是,在 MATLAB® 中,多通道图像排列为宽度×高度×通道数数组。要重构数据以使通道处于第三个维度中,请使用辅助函数 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 个类之一。获取具有对应 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
imshow(B)
title('Training Labels')
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');

为训练创建随机补片提取数据存储

使用随机补片提取数据存储将训练数据馈送给网络。此数据存储从一个包含真实值图像的图像数据存储和包含像素标签的像素标签数据存储中提取多个对应的随机补片。补片是防止大图像耗尽内存和有效增加可用训练数据量的常用方法。

首先将来自 '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)。每个小批量包含 16 个大小为 256×256 像素的补片。在轮的每次迭代中提取一千个小批量。

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

随机补片提取数据存储 dsTrain 在一轮训练的每次迭代中向网络提供小批量数据。预览数据存储以探查数据。

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 进行了修改以在卷积中使用零填充,从而使卷积的输入和输出具有相同的大小。使用辅助函数 createUnet 创建一个具有几个预选超参数的 U-Net。此函数作为支持文件包含在本示例中。

inputTileSize = [256,256,6];
lgraph = createUnet(inputTileSize);
disp(lgraph.Layers)
  58×1 Layer array with layers:

     1   'ImageInputLayer'                        Image Input                  256×256×6 images with 'zerocenter' normalization
     2   'Encoder-Section-1-Conv-1'               Convolution                  64 3×3×6 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 3×3×64 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                  2×2 max pooling with stride [2  2] and padding [0  0  0  0]
     7   'Encoder-Section-2-Conv-1'               Convolution                  128 3×3×64 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 3×3×128 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                  2×2 max pooling with stride [2  2] and padding [0  0  0  0]
    12   'Encoder-Section-3-Conv-1'               Convolution                  256 3×3×128 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 3×3×256 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                  2×2 max pooling with stride [2  2] and padding [0  0  0  0]
    17   'Encoder-Section-4-Conv-1'               Convolution                  512 3×3×256 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 3×3×512 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                  2×2 max pooling with stride [2  2] and padding [0  0  0  0]
    23   'Mid-Conv-1'                             Convolution                  1024 3×3×512 convolutions with stride [1  1] and padding [1  1  1  1]
    24   'Mid-ReLU-1'                             ReLU                         ReLU
    25   'Mid-Conv-2'                             Convolution                  1024 3×3×1024 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 2×2×1024 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 3×3×1024 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 3×3×512 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 2×2×512 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 3×3×512 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 3×3×256 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 2×2×256 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 3×3×256 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 3×3×128 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 2×2×128 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 3×3×128 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 3×3×64 convolutions with stride [1  1] and padding [1  1  1  1]
    55   'Decoder-Section-4-ReLU-2'               ReLU                         ReLU
    56   'Final-ConvolutionLayer'                 Convolution                  18 1×1×64 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) 优化来训练网络。使用 trainingOptions 函数指定 SGDM 的超参数设置。

训练深度网络是很费时间的。通过指定高学习率可加快训练速度。然而,这可能会导致网络的梯度爆炸或不受控制地增长,阻碍网络训练成功。要将梯度保持在有意义的范围内,请通过将 'GradientThreshold' 指定为 0.05 来启用梯度裁剪,并指定 'GradientThresholdMethod' 使用梯度的 L2-范数。

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

训练网络

默认情况下,示例使用 downloadTrainedUnet 辅助函数为该数据集下载 U-Net 的预训练版本。此函数作为支持文件包含在本示例中。借助预训练网络,您无需等待训练完成即可运行整个示例。

要训练网络,请将以下代码中的 doTraining 变量设置为 true。使用 trainNetwork 函数训练模型。

在 GPU 上(如果有)进行训练。使用 GPU 需要 Parallel Computing Toolbox™ 和支持 CUDA® 的 NVIDIA® GPU。有关详细信息,请参阅GPU Support by Release (Parallel Computing Toolbox)。在 NVIDIA Titan X 上训练大约需要 20 个小时。

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 对多光谱图像进行语义分割。

对测试数据预测结果

要基于经过训练的网络执行前向传导,请对验证数据集使用辅助函数 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%。

计算植被覆盖率

此示例的最终目标是计算多光谱图像中的植被覆盖率。

求已加植被标签的像素的数量。标签 ID 2 ("Trees")、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."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.

另请参阅

| | (Image Processing Toolbox) | (Computer Vision Toolbox) | (Computer Vision Toolbox) | (Computer Vision Toolbox) | | (Image Processing Toolbox) | (Computer Vision Toolbox)

相关主题

外部网站