航空 LiDAR データの地形分類

この例では、航空 LiDAR データの地形を地面、建物、および植生としてセグメント化して分類する方法を示します。この例では航空 LiDAR システムで取得された LAZ ファイルを入力として使用します。最初に、LAZ ファイルの点群データを地面と地面以外の点に分類します。その後、地面以外の点を法線と曲率の特徴に基づいて建物と植生の点にさらに分類します。次の図にプロセスの概要を示します。

データの読み込みと可視化

Open Topography Dataset [1] から取得した LAZ ファイル aerialLidarData.laz から、点群データと対応するグラウンド トゥルース ラベルを読み込みます。この点群は、地面、建物、および植生を含むさまざまなクラスで構成されています。lasFileReaderオブジェクトのオブジェクト関数readPointCloudを使用して、点群データと対応するグラウンド トゥルース ラベルをワークスペースに読み込みます。関数pcshowを使用して、グラウンド トゥルース ラベルに従って色分けされた点群を可視化します。

lazfile = fullfile(toolboxdir('lidar'),'lidardata','las','aerialLidarData.laz');
% Read LAZ data from file
lazReader = lasFileReader(lazfile);
% Read point cloud and corresponding ground truth labels
[ptCloud,pointAttributes] = readPointCloud(lazReader, ...
    'Attributes','Classification');
grdTruthLabels = pointAttributes.Classification;
% Visualize the input point cloud with corresponding ground truth labels
figure
pcshow(ptCloud.Location,grdTruthLabels)
title('Aerial Lidar Data with Ground Truth')

地面の分類

"地面の分類" は、入力点群を地面と地面以外にセグメント化する前処理の手順です。関数segmentGroundSMRFを使用して、LAZ ファイルから読み込んだデータを地面と地面以外の点にセグメント化します。

[groundPtsIdx,nonGroundPtCloud,groundPtCloud] = segmentGroundSMRF(ptCloud);
% Visualize ground and non-ground points in green and magenta, respectively
figure
pcshowpair(nonGroundPtCloud,groundPtCloud)
title('Classified Ground and Non-Ground Points')

特徴抽出

関数 helperExtractFeatures を使用して、点群から特徴を抽出します。この例には、すべての補助関数がサポート ファイルとして添付されています。この補助関数は、点群に含まれる各点の法線と曲率の値を推定します。これらの特徴は、各点を近傍の点と関連付けて各点の基となる構造の情報を提供します。

近傍と見なす数を指定できます。近傍の数が小さすぎると、補助関数で植生の点が過剰にクラスター化されます。近傍の数が大きすぎると、建物と植生の間の明確な境界がなくなり、建物の点の近くにある植生の点が誤って分類されます。

neighbors = 10;  
[normals,curvatures,neighInds] = helperExtractFeatures(nonGroundPtCloud, ...
    neighbors);

建物と植生の分類

補助関数では、法線と曲率の変化を使用して建物と植生を区別します。建物は植生に比べて平面的であるため、建物に属する点は曲率の変化や近傍間の法線の相対差が小さくなります。植生の点は、建物に比べて散在しているため、曲率の変化が大きくなります。関数 helperClassify は、地面以外の点を建物と植生に分類します。この補助関数は、次の基準に基づいて、該当する点を建物として分類します。

• 各点の曲率が小さく、curveThresh で指定された曲率のしきい値以内でなければならない。

• 近傍点の法線が類似していなければならない。近傍の法線間の余弦類似性が normalThresh で指定された法線のしきい値より大きくなければならない。

上記の基準を満たさない点は植生としてマークされます。この補助関数は、植生に属する点を 1、建物に属する点を 2 としてラベル付けします。

% Specify the normal threshold and curvature threshold
normalThresh = 0.85;
curveThresh = 0.02;
% Classify the points into building and vegetation
labels = helperClassify(normals,curvatures,neighInds, ...
    normalThresh,curveThresh);

建物と植生のクラス ラベルをグラウンド トゥルース ラベル データから抽出します。LAZ ファイルには多くのクラスが含まれているため、最初に地面、建物、および植生のクラスを分離する必要があります。分類ラベルは ASPRS の LAZ ファイル形式の標準に準拠しています。

• 分類値 2 — 地面の点を表す

• 分類値 3、4、5 — 低植生、中植生、高植生の点を表す

• 分類値 6 — 建物の点を表す

maskData を定義して、地面、建物、および植生に属する点を入力点群から抽出します。

maskData = grdTruthLabels>=2 & grdTruthLabels<=6;

grdTruthLabels として指定された入力点群のグラウンド トゥルース ラベルを変更します。

% Compress low, medium, and high vegetation to a single value
grdTruthLabels(grdTruthLabels>=3 & grdTruthLabels<=5) = 4;
% Update grdTruthLabels for metrics calculation
grdTruthLabels(grdTruthLabels == 2) = 1;    
grdTruthLabels(grdTruthLabels == 4) = 2;
grdTruthLabels(grdTruthLabels == 6) = 3;

前の分類の手順から取得した予測ラベルを estimatedLabels に格納します。

estimatedLabels = zeros(ptCloud.Count,1);
estimatedLabels(groundPtsIdx) = 1; 
estimatedLabels(labels == 1) = 2;
estimatedLabels(labels == 2) = 3;

地面、建物、および植生に属するラベルを抽出します。

grdTruthLabels = grdTruthLabels(maskData);
estimatedLabels = estimatedLabels(maskData);

グラウンド トゥルース ラベルと推定ラベルで地形を可視化します。

ptCloud = select(ptCloud,maskData);
hFig = figure('Position',[0 0 900 400]);
axMap1 = subplot(1,2,1,'Color','black','Parent',hFig);
axMap1.Position = [0 0.2 0.5 0.55];
pcshow(ptCloud.Location,grdTruthLabels,'Parent',axMap1)
axis off
title(axMap1,'Aerial Lidar Data with Ground Truth Labels')
axMap2 = subplot(1,2,2,'Color','black','Parent',hFig);
axMap2.Position = [0.5,0.2,0.5,0.55];
pcshow(ptCloud.Location,estimatedLabels,'Parent',axMap2)
axis off
title(axMap2,'Aerial Lidar Data with Classified Labels')

検証

与えられた点群の全体的な精度とクラスの精度、Intersection-over-Union (IoU)、および加重 IoU を計算して分類を検証します。

confusionMatrix = segmentationConfusionMatrix(estimatedLabels,double(grdTruthLabels));
ssm = evaluateSemanticSegmentation({confusionMatrix}, ...
    {'Ground' 'Vegetation' 'Building'},'Verbose',0);
disp(ssm.DataSetMetrics)

    GlobalAccuracy    MeanAccuracy    MeanIoU    WeightedIoU
    ______________    ____________    _______    ___________

       0.71112          0.54571       0.32094      0.66208  

disp(ssm.ClassMetrics)

                  Accuracy      IoU   
                  ________    ________

    Ground        0.72222      0.72222
    Vegetation    0.80547      0.20532
    Building      0.10944     0.035268

参考

関数

readPointCloud | segmentGroundSMRF| pcnormals | pcshow | pcshowpair | segmentationConfusionMatrix | evaluateSemanticSegmentation

オブジェクト

lasFileReader

参考文献

[1] Starr, Scott. "Tuscaloosa, AL: Seasonal Inundation Dynamics and Invertebrate Communities." National Center for Airborne Laser Mapping, December 1, 2011. OpenTopography (https://doi.org/10.5069/G9SF2T3K)