Main Content

深層学習を使用したマルチスペクトル イメージのセマンティック セグメンテーション

この例では、U-Net を使用して 7 チャネルのマルチスペクトル イメージのセマンティック セグメンテーションを実行する方法を示します。

セマンティック セグメンテーションでは、イメージの各ピクセルにクラスでラベル付けします。セマンティック セグメンテーションの適用例の 1 つに、森林被覆の時間的変化である森林伐採の追跡があります。環境保護機関は、森林伐採を追跡し、地域の環境的生態学的健全性を評価し測定します。

深層学習ベースのセマンティック セグメンテーションにより、高解像度の航空写真から植被を正確に測定できます。課題の 1 つに、視覚的に類似した特性のクラスを切り分けること、たとえば緑のピクセルを草、低木または樹木として区別することがあります。分類の精度を高めるために、一部のデータ セットには各ピクセルに関する追加情報を提供するマルチスペクトル イメージが含まれています。たとえば、ハームリン ビーチ州立公園のデータ セットでは、クラスをより明確に分離する 3 つの近赤外チャネルでカラー イメージが補完されています。

この例では、まず、事前学習済みの U-Net を使用してセマンティック セグメンテーションを実行する方法を示し、次に、セグメンテーションの結果を使用して植被率を計算する方法を示します。さらに、オプションとして、パッチベースの学習手法を使用して、Hamlin Beach State Park データ セットで U-Net ネットワークに学習させることもできます。

事前学習済みの U-Net のダウンロード

学習済みネットワークとデータ セットの目的の場所として dataDir を指定します。

dataDir = fullfile(tempdir,"rit18_data");

事前学習済みの U-Net ネットワークをダウンロードします。

trainedNet_url = "https://www.mathworks.com/supportfiles/"+ ...
    "vision/data/trainedMultispectralUnetModel.zip";
downloadTrainedNetwork(trainedNet_url,dataDir);
load(fullfile(dataDir,"trainedMultispectralUnetModel", ...
    "trainedMultispectralUnetModel.mat"));

データセットのダウンロード

この例では、ネットワークに学習させるために、高解像度のマルチスペクトル データ セットを使用します [1]。このイメージ セットはニューヨーク州のハームリン ビーチ州立公園の上空でドローンを使用して撮影されました。このデータには、18 個のオブジェクト クラス ラベルの付いた、ラベル付き学習セット、検証セットおよびテスト セットが含まれます。データ ファイルのサイズは 3.0 GB です。

補助関数 downloadHamlinBeachMSIData を使用して、データ セットの MAT ファイル バージョンをダウンロードします。この関数は、この例にサポート ファイルとして添付されています。

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              

マルチスペクトル イメージ データは numChannels x width x height 配列に配置されます。ただし、MATLAB® では、マルチチャネル イメージは width x height x numChannels 配列に配置されます。チャネルが 3 番目の次元になるようにデータを形状変更するには、関数 permute を使用します。

train_data = permute(train_data,[2 3 1]);
val_data = permute(val_data,[2 3 1]);
test_data = permute(test_data,[2 3 1]);

データが正しい構造体であることを確認します。

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              

マルチスペクトル データの可視化

各スペクトル バンドの中心をナノメートル単位で表示します。

disp(band_centers)
   490   550   680   720   800   900

このデータ セットでは、RGB カラー チャネルがそれぞれ 3 番目、2 番目、および 1 番目のイメージ チャネルになります。学習イメージ、検証イメージおよびテスト イメージの RGB 成分をモンタージュとして表示します。イメージを画面上で明るく表示するには、関数histeqを使用してヒストグラム均等化を行います。

rgbTrain = histeq(train_data(:,:,[3 2 1]));
rgbVal = histeq(val_data(:,:,[3 2 1]));
rgbTest = histeq(test_data(:,:,[3 2 1]));

montage({rgbTrain,rgbVal,rgbTest},BorderSize=10,BackgroundColor="white")
title("RGB Component of Training, Validation, and Test Image (Left to Right)")

データの 4 番目、5 番目、および 6 番目のチャネルは近赤外帯域に対応します。学習イメージのこれら 3 つのチャネルのヒストグラムを均等化し、これらのチャネルをモンタージュとして表示します。チャネルは、熱の痕跡に基づいてイメージの異なる成分を強調します。たとえば、4 番目のチャネルでは、他の 2 つの赤外線チャネルよりも木が暗くなります。

ir4Train = histeq(train_data(:,:,4));
ir5Train = histeq(train_data(:,:,5));
ir6Train = histeq(train_data(:,:,6));

montage({ir4Train,ir5Train,ir6Train},BorderSize=10,BackgroundColor="white")
title("Infrared Channels 4, 5, and 6 (Left to Right) of Training Image ")

データの 7 番目のチャネルは、有効なセグメンテーション領域を示すバイナリ マスクです。学習イメージ、検証イメージおよびテスト イメージのマスクを表示します。

maskTrain = train_data(:,:,7);
maskVal = val_data(:,:,7);
maskTest = test_data(:,:,7);

montage({maskTrain,maskVal,maskTest},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)

この例は、イメージを植生と非植生の 2 つのクラスにセグメント化することを目的としています。ターゲット クラス名を定義します。

classNames = ["NotVegetation" "Vegetation"];

18 個の元のクラスを、学習データ用と検証データ用に 2 つのターゲット クラスにグループ化します。"Vegetation" は、クラス ID 2、13、14 をもつ元のクラス "Tree"、"Other Vegetable"、および "Grass" の組み合わせです。クラス ID 0 の元のクラス "Other Class/Image Border" は背景クラスに属します。他のすべての元のクラスは、ターゲット ラベル "NotVegetation" に属します。

vegetationClassIDs = [2 13 14];
nonvegetationClassIDs = setdiff(1:length(classes),vegetationClassIDs);

labelsTrain = zeros(size(train_labels),"uint8");
labelsTrain(ismember(train_labels,nonvegetationClassIDs)) = 1;
labelsTrain(ismember(train_labels,vegetationClassIDs)) = 2;

labelsVal = zeros(size(val_labels),"uint8");
labelsVal(ismember(val_labels,nonvegetationClassIDs)) = 1;
labelsVal(ismember(val_labels,vegetationClassIDs)) = 2;

グラウンド トゥルース検証ラベルを PNG ファイルとして保存します。この例では、このファイルを使用して精度メトリクスを計算します。

imwrite(labelsVal,"gtruth.png");

ヒストグラム均等化を行った RGB 学習イメージにラベルを重ね合わせます。カラー バーをイメージに追加します。

cmap = [1 0 1;0 1 0];
B = labeloverlay(rgbTrain,labelsTrain,Transparency=0.8,Colormap=cmap);
imshow(B,cmap)
title("Training Labels")
numClasses = numel(classNames);
ticks = 1/(numClasses*2):1/numClasses:1;
colorbar(TickLabels=cellstr(classNames),Ticks=ticks,TickLength=0,TickLabelInterpreter="none");

テスト イメージでのセマンティック セグメンテーションの実行

このイメージのサイズでは、イメージ全体を一度にセグメント化することができません。代わりに、ブロック化イメージの手法を使用してイメージをセグメント化します。この手法では、一度に 1 つのデータ ブロックを読み込んで処理するため、非常に大きなファイルにも対応できます。

関数blockedImageを使用して、テスト データの 6 つのスペクトル チャネルを含むブロック化イメージを作成します。

patchSize = [1024 1024];
bimTest = blockedImage(test_data(:,:,1:6),BlockSize=patchSize);

関数semanticseg (Computer Vision Toolbox)を使用して、データ ブロックをセグメント化します。関数 apply を使用して、ブロック化イメージ内のすべてのブロックに対して関数 sematicseg を呼び出します。

bimSeg = apply(bimTest,@(bs)semanticseg(bs.Data,net,Outputtype="uint8"),...
    PadPartialBlocks=true,PadMethod=0);

関数 gather を使用して、セグメント化されたすべてのブロックを 1 つのイメージとしてワークスペースにまとめます。

segmentedImage = gather(bimSeg);

セグメンテーションの有効な部分だけを抽出するには、セグメント化されたイメージに検証データのマスク チャネルを乗算します。

segmentedImage = segmentedImage .* uint8(maskTest~=0);
imshow(segmentedImage,[])
title("Segmented Image")

セマンティック セグメンテーションの出力にはノイズが多く含まれます。イメージの後処理を実行して、ノイズと散在ピクセルを取り除きます。関数medfilt2を使用して、セグメンテーションからごま塩ノイズを取り除きます。ノイズが除去され、セグメント化されたイメージを表示します。

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

ヒストグラム均等化を行った RGB 検証イメージにセグメント化されたイメージを重ね合わせます。

B = labeloverlay(rgbTest,segmentedImage,Transparency=0.8,Colormap=cmap);
imshow(B,cmap)
title("Labeled Segmented Image")
colorbar(TickLabels=cellstr(classNames),Ticks=ticks,TickLength=0,TickLabelInterpreter="none");

植被率の計算

セマンティック セグメンテーションの結果を使用して、関連する生態系に関する質問に答えることができます。たとえば、植生に覆われている土地の割合を知りたいとします。この質問に答えるには、セグメント化されたテスト イメージ内で植生としてラベル付けされたピクセルの数を求めます。また、セグメント化されたイメージ内の非ゼロのピクセル数を数えることにより、ROI 内のピクセルの合計数を求めます。

vegetationPixels = ismember(segmentedImage(:),vegetationClassIDs);
numVegetationPixels = sum(vegetationPixels(:));
numROIPixels = nnz(segmentedImage);

植生のピクセルの数を ROI 内のピクセルの数で除算して、植被率を計算します。

percentVegetationCover = (numVegetationPixels/numROIPixels)*100;
disp("The percentage of vegetation cover is "+percentVegetationCover+"%");
The percentage of vegetation cover is 65.8067%

この例の残りの部分では、Hamlin Beach データ セットで U-Net に学習させる方法を示します。

学習用のブロック化イメージ データストアの作成

ブロック化イメージ データストアを使用して、ネットワークに学習データを供給します。このデータストアは、グラウンド トゥルース イメージとピクセル ラベル データを含むイメージ データストアとピクセル ラベル データストアから、対応する複数のパッチを抽出します。

学習イメージ、学習ラベル、マスクをブロック化されたイメージとして読み取ります。

inputTileSize = [256 256];
bim = blockedImage(train_data(:,:,1:6),BlockSize=inputTileSize);
bLabels = blockedImage(labelsTrain,BlockSize=inputTileSize);
bmask = blockedImage(maskTrain,BlockSize=inputTileSize);

マスクとオーバーラップするイメージ データのブロックを選択します。

overlapPct = 0.185;
blockOffsets = round(inputTileSize.*overlapPct);
bls = selectBlockLocations(bLabels, ...
    BlockSize=inputTileSize,BlockOffsets=blockOffsets, ...
    Masks=bmask,InclusionThreshold=0.95);

ラベルを one-hot 符号化します。

labelsTrain1hot = onehotencode(labelsTrain,3,ClassNames=1:2);
labelsTrain1hot(isnan(labelsTrain1hot)) = 0;
bLabels = blockedImage(labelsTrain1hot,BlockSize=inputTileSize);

関数blockedImageDatastoreを使用して、ブロック化イメージ データストアにデータを書き込みます。

bimds = blockedImageDatastore(bim,BlockLocationSet=bls,PadMethod=0);
bimdsLabels = blockedImageDatastore(bLabels,BlockLocationSet=bls,PadMethod=0);

2 つのブロック化イメージ データストアからCombinedDatastoreを作成します。

dsTrain = combine(bimds,bimdsLabels);

ブロック化イメージ データストア dsTrain は、エポックの各反復でデータのミニバッチをネットワークに渡します。データストアをプレビューしてデータを調査します。

preview(dsTrain)
ans=1×2 cell array
    {256×256×6 uint16}    {256×256×2 double}

U-Net ネットワーク層の作成

この例では U-Net ネットワークのバリエーションを使用します。U-Net では、最初の一連の畳み込み層に最大プーリング層が点在し、入力イメージの解像度を逐次下げていきます。これらの層に、一連の畳み込み層が続き、その中にアップサンプリング演算処理が点在し、入力イメージの解像度を逐次上げていきます [2]。U-Net の名前は、このネットワークが文字「U」のように対称の形状で描けることに由来しています。

U-Net のハイパーパラメーターを指定します。入力深度はハイパースペクトル チャネル数 6 です。

inputDepth = 6;
encoderDepth = 4;
convFilterSize = 3;
upconvFilterSize = 2;

関数 blockedNetwork を使用して、層の複数の反復ブロックで構成される符号化器モジュールを作成します。補助関数 encoderBlockMultispectralUNet は、符号化器の層からなる単一のブロックを作成します。これは、この例にサポート ファイルとして添付されています。

encoderBlockFcn = @(block) ...
    encoderBlockMultispectralUNet(block,inputDepth,convFilterSize,encoderDepth);
encoder = blockedNetwork(encoderBlockFcn,encoderDepth,NamePrefix="encoder_");

関数 blockedNetwork を使用して、層の複数の反復ブロックで構成される復号化器モジュールを作成します。補助関数 decoderBlockMultispectralUNet は、復号化器の層からなる単一のブロックを作成します。これは、この例にサポート ファイルとして添付されています。

decoderBlockFcn = @(block) ...
    decoderBlockMultispectralUNet(block,convFilterSize,upconvFilterSize);
decoder = blockedNetwork(decoderBlockFcn,encoderDepth,NamePrefix="decoder_");

補助関数 bridgeBlockMultispectralUNet を使用してブリッジ層を定義します。これは、この例にサポート ファイルとして添付されています。

bridge = bridgeBlockMultispectralUNet(convFilterSize,encoderDepth);

出力層を定義します。

final = [
    convolution2dLayer(1,numClasses,Padding="same")
    softmaxLayer];

関数 encoderDecoderNetwork を使用して、符号化器モジュール、ブリッジ、復号化器モジュール、および最終層を接続します。スキップ接続を追加します。

skipConnectionNames = [
    "encoder_Block1Layer5","decoder_Block4Layer2";
    "encoder_Block2Layer5","decoder_Block3Layer2";
    "encoder_Block3Layer5","decoder_Block2Layer2";
    "encoder_Block4Layer5","decoder_Block1Layer2"];
unet = encoderDecoderNetwork([inputTileSize inputDepth],encoder,decoder, ...
    OutputChannels=numClasses, ...
    SkipConnectionNames=skipConnectionNames, ...
    SkipConnections="concatenate", ...
    LatentNetwork=bridge, ...
    FinalNetwork=final);

学習オプションの選択

モーメンタム項付き確率的勾配降下 (SGDM) 最適化を使用してネットワークに学習させます。関数trainingOptions (Deep Learning Toolbox)を使用して SGDM 用ハイパーパラメーター設定を指定します。勾配のクリップを有効にするために、名前と値の引数 GradientThreshold0.05 に指定し、GradientThresholdMethod を指定して勾配の L2 ノルムを使用します。

maxEpochs = 150;
minibatchSize = 16;

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

ネットワークの学習

ネットワークに学習させるには、次のコードで変数 doTrainingtrue に設定します。関数trainnet (Deep Learning Toolbox)を使用して、モデルに学習させます。マスクされていないピクセルのみのクロス エントロピー損失を計算するカスタム損失関数 modelLoss を指定します。このカスタム損失関数は、この例の終わりで定義されています。

GPU が利用できる場合、GPU で学習を行います。GPU を使用するには、Parallel Computing Toolbox™、および CUDA® 対応の NVIDIA® GPU が必要です。詳細については、GPU 計算の要件 (Parallel Computing Toolbox)を参照してください。

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

セグメンテーションの精度の評価

検証データをセグメント化します。

関数blockedImageを使用して、検証データの 6 つのスペクトル チャネルを含むブロック化イメージを作成します。

bimVal = blockedImage(val_data(:,:,1:6),BlockSize=patchSize);

関数semanticseg (Computer Vision Toolbox)を使用して、データ ブロックをセグメント化します。関数 apply を使用して、ブロック化イメージ内のすべてのブロックに対して関数 sematicseg を呼び出します。

bimSeg = apply(bimVal,@(bs)semanticseg(bs.Data,net,Outputtype="uint8"),...
    PadPartialBlocks=true,PadMethod=0);

関数 gather を使用して、セグメント化されたすべてのブロックを 1 つのイメージとしてワークスペースにまとめます。

segmentedImage = gather(bimSeg);

セグメント化されたイメージを PNG ファイルとして保存します。

imwrite(segmentedImage,"results.png");

関数pixelLabelDatastore (Computer Vision Toolbox)を使用して、セグメンテーションの結果とグラウンド トゥルース ラベルを読み込みます。

pixelLabelIDs = [1 2];
pxdsResults = pixelLabelDatastore("results.png",classNames,pixelLabelIDs);
pxdsTruth = pixelLabelDatastore("gtruth.png",classNames,pixelLabelIDs);

関数evaluateSemanticSegmentation (Computer Vision Toolbox)を使用して、セマンティック セグメンテーションの精度を測定します。グローバル精度のスコアは、正しく分類されたピクセルの割合が 96% を超えていることを示しています。

ssm = evaluateSemanticSegmentation(pxdsResults,pxdsTruth);
Evaluating semantic segmentation results
----------------------------------------
* Selected metrics: global accuracy, class accuracy, IoU, weighted IoU, BF score.
* Processed 1 images.
* Finalizing... Done.
* Data set metrics:

    GlobalAccuracy    MeanAccuracy    MeanIoU    WeightedIoU    MeanBFScore
    ______________    ____________    _______    ___________    ___________

       0.96875          0.96762       0.93914      0.93931        0.79113  

補助関数

関数 modelLoss は、イメージのマスクされていないすべてのピクセルに対してクロス エントロピー損失を計算します。

function loss = modelLoss(y,targets)
    mask = ~isnan(targets);
    targets(isnan(targets)) = 0;
    loss = crossentropy(y,targets,Mask=mask);
end

参考文献

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

参考

| | (Deep Learning Toolbox) | (Deep Learning Toolbox) | | | | | (Computer Vision Toolbox) | (Computer Vision Toolbox) | (Computer Vision Toolbox)

関連するトピック

外部の Web サイト