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

此示例使用：

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

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

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

此示例首先说明如何使用预训练的 U-Net 执行语义分割，然后使用分割结果来计算植被覆盖率。然后，您可以选择使用基于补片的训练方法基于 Hamlin Beach State Parck 数据集训练 U-Net 网络。

下载数据集

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

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

dataDir = fullfile(tempdir,"rit18_data"); 
downloadHamlinBeachMSIData(dataDir);

加载该数据集。

load(fullfile(dataDir,"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

将训练数据保存为一个 MAT 文件，将训练标签保存为一个 PNG 文件。这有助于在训练期间使用 imageDatastore 和 pixelLabelDatastore 加载训练数据。

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

可视化多光谱数据

在此数据集中，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, Validation, and Test Image (Left to Right)")

将训练数据的后三个经过直方图均衡化处理的通道显示为蒙太奇形式。这些通道对应于近红外波段，并基于其热特征突出显示图像的不同分量。例如，第二个通道图像中心附近的树比其他两个通道中的树显示的细节要多。

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

通道 7 是指示有效分割区域的掩膜。显示训练、验证和测试图像的掩膜。

figure
montage(...
    {train_data(:,:,7),val_data(:,:,7),test_data(:,:,7)}, ...
    BorderSize=10,BackgroundColor="white")
title("Mask of Training, Validation, and Test Image (Left to 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)

执行语义分割

下载预训练的 U-Net 网络。

trainedUnet_url = "https://www.mathworks.com/supportfiles/vision/data/multispectralUnet.mat";
downloadTrainedNetwork(trainedUnet_url,dataDir);
load(fullfile(dataDir,"multispectralUnet.mat"));

为了在经过训练的网络上执行语义分割，请对验证数据使用 segmentMultispectralImage 辅助函数。此函数作为支持文件包含在本示例中。segmentMultispectralImage 函数使用 semanticseg (Computer Vision Toolbox) 函数对图像补片执行分割。需要处理补片，因为图像的大小导致无法一次性处理整个图像。

predictPatchSize = [1024 1024];
segmentedImage = segmentMultispectralImage(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 Segmented Image")
colorbar(TickLabels=cellstr(classNames),Ticks=ticks,TickLength=0,TickLabelInterpreter="none");
colormap(cmap)

计算植被覆盖率

语义分割结果可用于解决相关生态学问题。例如，植被覆盖的陆地面积百分比是多少？要回答此问题，请找出标注为植被的像素数。标签 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%.

该示例的其余部分说明如何基于 Hamlin Beach 数据集训练 U-Net。

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

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

首先从 imageDatastore 中的 "train_data.mat" 加载训练图像。由于 MAT 文件格式是非标准图像格式，您必须使用 MAT 文件读取器来读取图像数据。您可以使用 MAT 文件读取器辅助函数 matRead6Channels，它从训练数据中提取前六个通道，省略包含掩膜的最后一个通道。此函数作为支持文件包含在本示例中。

imds = imageDatastore("train_data.mat",FileExtensions=".mat",ReadFcn=@matRead6Channels);

创建一个 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'               2-D 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'               2-D 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'              2-D Max Pooling              2×2 max pooling with stride [2  2] and padding [0  0  0  0]
     7   'Encoder-Section-2-Conv-1'               2-D 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'               2-D 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'              2-D Max Pooling              2×2 max pooling with stride [2  2] and padding [0  0  0  0]
    12   'Encoder-Section-3-Conv-1'               2-D 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'               2-D 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'              2-D Max Pooling              2×2 max pooling with stride [2  2] and padding [0  0  0  0]
    17   'Encoder-Section-4-Conv-1'               2-D 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'               2-D 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'              2-D Max Pooling              2×2 max pooling with stride [2  2] and padding [0  0  0  0]
    23   'Mid-Conv-1'                             2-D 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'                             2-D 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'               2-D 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'               2-D 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'               2-D 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'               2-D 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'               2-D 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'               2-D 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'               2-D 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'               2-D 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'               2-D 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'               2-D 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'               2-D 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'               2-D 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'                 2-D 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);

训练网络或下载预训练网络

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

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

doTraining = false; 
if doTraining
    net = trainNetwork(dsTrain,lgraph,options);
    modelDateTime = string(datetime("now",Format="yyyy-MM-dd-HH-mm-ss"));
    save(fullfile(dataDir,"multispectralUnet-"+modelDateTime+".mat"),"net");
end

评估分割准确度

分割验证数据。

segmentedImage = segmentMultispectralImage(val_data,net,predictPatchSize);

将分割后的图像和真实值标签另存为 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.90411

全局准确度分数表明，正确分类的像素稍大于 90%。

参考资料

[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.

[3] Kemker, Ronald, Carl Salvaggio, and Christopher Kanan."Algorithms for Semantic Segmentation of Multispectral Remote Sensing Imagery Using Deep Learning."ISPRS Journal of Photogrammetry and Remote Sensing, Deep Learning RS Data, 145 (November 1, 2018):60-77. https://doi.org/10.1016/j.isprsjprs.2018.04.014.