Семантическая сегментация Aerial Lidar с использованием глубокого обучения PointNet++

Этот пример использует:

В этом примере показано, как обучить нейронную сеть для глубокого обучения PointNet++ для выполнения семантической сегментации на данных о воздушном лидаре.

Данные лидар, полученные из бортовых лазерных систем сканирования, все чаще используются в таких приложениях, как топографическое отображение, создание и городов, измерение биомассы и управление стихийными бедствиями. Извлечение содержательной информации из этих данных требует семантической сегментации, процесса, в котором каждому пикселю в облаке точек присваивается метка класса.

В этом примере вы обучаете сеть PointNet++ для выполнения семантической сегментации с помощью набора данных [1] Dayton annotated lidar скан (DALES), который содержит сцены плотных, маркированных лидарными воздушными данными из городских, пригородных, сельских и коммерческих настроек. Из 40 сцен 29 сцен используются для обучения, а остальные 11 последовательностей используются для проверки. Набор данных обеспечивает семантические метки сегментации для 8, таких как создания, автомобили, грузовики, столбы, линии степеней, ограждения, земля и растительность.

Загрузка данных DALES

Набор данных DALES содержит 40 сцен данных о воздушном лидаре. Набор данных уже разделен на 29 сцен для обучения и 11 сцен для проверки. Каждый пиксель данных имеет метку класса. Следуйте инструкциям на веб-сайте DALES, чтобы загрузить набор данных в папку, заданную dataFolder переменная. Создайте папки для хранения обучающих и тестовых данных.

dataFolder = fullfile(tempdir,'DALES');
trainDataFolder = fullfile(dataFolder,'dales_las','train');
testDataFolder = fullfile(dataFolder,'dales_las','test');

Предварительный просмотр облака точек из обучающих данных.

lasReader = lasFileReader(fullfile(trainDataFolder,'5080_54435.las'));
[pc,attr] = readPointCloud(lasReader,"Attributes", "Classification");
labels = attr.Classification + 1;
classNames = [
    "background"
    "ground"
    "vegetation"
    "cars"
    "trucks"
    "powerlines"
    "poles"
    "fences"
    "buildings"
    ];
cmap = labelColorMap();
figure;
colormap(cmap);
ax = pcshow(pc.Location,labels);
labelColorbar(ax,cmap,classNames);
title("PointCloud with overlaid semantic labels");

Предварительная обработка данных

Каждое облако точек в наборе данных DALES занимает площадь 500 на 500 метров, что намного больше типичной площади, охватываемой наземным вращающимся лидаром облаков точек. Для эффективной обработки памяти разделите облако точек на небольшие, неперекрывающиеся сетки. Кроме того, сосредоточьте внимание на трех доминирующих классах - земле, созданиях и растительности - чтобы избежать проблем из-за классового дисбаланса.

Используйте helperCropPointCloudsAndMergeLabels функция, присоединенная к этому примеру как вспомогательному файлу, для:

Обрезать облака точек на неперекрывающиеся сетки размером 50 на 50 метров.
Объедините все классы, кроме доминирующих классов (землю, создания и растительность), в один фоновый класс, так чтобы данные содержали только четыре класса.
Сохраните обрезанные сетки и семантические метки как файлы PCD и PNG, соответственно.

Определите размерности сетки и установите фиксированное число точек на сетку, чтобы обеспечить более быстрое обучение.

gridSize = [50,50];
numPoints = 115000;

Если вы используете этот рабочий процесс для собственных данных, задайте writeFiles «ложь», если обучающие данные разделены на сетки. Обратите внимание, что обучающие данные должны быть в формате, поддерживаемом pcread функция.

writeFiles = true;
[pcCropTrainPath,labelsCropTrainPath] = helperCropPointCloudsAndMergeLabels(gridSize,trainDataFolder,numPoints,writeFiles);
[pcCropTestPath,labelsCropTestPath] = helperCropPointCloudsAndMergeLabels(gridSize,testDataFolder,numPoints,writeFiles);

Обработка может занять некоторое время. Код приостанавливает выполнение MATLAB ® до завершения обработки.

Создайте объекты Datastore для обучения

Создайте fileDatastore объект для загрузки файлов PCD с помощью pcread функция.

ldsTrain = fileDatastore(pcCropTrainPath,'ReadFcn',@(x) pcread(x));

Использование pixelLabelDatastore объект для хранения пиксельных меток из изображений меток пикселей. Объект сопоставляет каждую пиксельную метку с именем класса. В этом примере земля, растительность и создания являются единственными интересующими классами; все другие пиксели являются фоном. Задайте эти классы и присвойте уникальный идентификатор метки каждому классу.

classNames = classNames([1,2,3,9],:)

classNames = 4×1 string
    "background"
    "ground"
    "vegetation"
    "buildings"

numClasses = numel(classNames);

% Specify label IDs from 1 to the number of classes.
labelIDs = 1 : numClasses;
pxdsTrain = pixelLabelDatastore(labelsCropTrainPath,classNames,labelIDs);

Загрузка и отображение облаков точек.

ptcld = preview(ldsTrain);
labels = preview(pxdsTrain);
cmap = cmap([1,2,3,9],:);
figure;
ax1 = pcshow(ptcld.Location,uint8(labels));
labelColorbar(ax1,cmap,classNames);
title("Cropped point cloud with overlaid semantic labels");

Используйте normalizePointCloud функция, заданная в конце примера, чтобы нормализовать облако точек между 0 и 1.

ldsTransformed = transform(ldsTrain,@(x) normalizePointCloud(x));

Используйте onehotEncodeLabels функция, заданная в конце примера, для преобразования меток в представление с одним горячим кодированием для настройки сети сегментации с помощью пользовательского цикла обучения.

pxdsTransformed = transform(pxdsTrain,@(x) onehotEncodeLabels(x,classNames));

Используйте combine функция для объединения облаков точек и меток пикселей в одном datastore для обучения.

dsTrain = combine(ldsTransformed,pxdsTransformed);

Задайте модель PointNet++

Модель сегментации PointNet++ [2] состоит из двух основных компонентов:

Установите абстракционные модули
Модули распространения функций

Серия модулей абстракции набора постепенно подсеменяет интересующие точки иерархически группированием точек и использует пользовательскую архитектуру PointNet, чтобы закодировать точки в векторы функций. Поскольку семантические задачи сегментации требуют функций точек для всех исходных точек, этот пример использует серию модулей распространения функций для иерархической интерполяции функций к исходным точкам с помощью схемы интерполяции на основе обратного расстояния.

Модуль абстракции набора реализован с помощью sampleAndGroup и sharedMLP вспомогательные функции и модуль распространения признаков реализован с помощью featurePropagation и sharedMLP вспомогательные функции. The sharedMLP модуль является единственным обучаемым модулем, который состоит из слоев свертки и нормализации образца. Используйте вспомогательную функцию pointNetplusModel для определения всей модели сегментации PointNet++. Вспомогательные функции перечислены в конце этого примера.

Задайте параметры модели PointNet++

The sampleAndGroup функция модуля абстракции набора параметризована количеством кластеров, количеством выборок на кластер и радиусом каждого кластера, тогда как sharedMLP функция реализует пользовательскую архитектуру PointNet, которая параметризована с помощью количества входа каналов и скрытых размеров каналов.

В этом примере используются значения гиперзначений параметров, настроенные на наборе данных DALES. Если вы хотите применить PointNet++ к другому набору данных, необходимо выполнить дополнительную настройку гиперпараметра.

Установите параметры первого модуля абстракции набора. Установите количество кластеров равное 1024, количество выборок в каждом кластере равное 128 и радиус каждого кластера равный 0,1. Также установите размер входного канала модели PointNet равным трем, а размер скрытого канала равным 32, 32 и 64.

inputChannelSize = 3;
hiddenChannelSize = [32,32,64];
nontrainable.nClusters1 = 1024;
nontrainable.nSamples1 = 128;
nontrainable.radius1 = 0.1;
trainable.SharedMLP1 = initializeSharedMLP(inputChannelSize,hiddenChannelSize);

Установите параметры второго модуля абстракции набора. Установите количество кластеров равное 256, количество выборок в каждом кластере равное 64 и радиус каждого кластера равный 0,2. Также установите размер входного канала модели PointNet равным 67, а размер скрытого канала равным 64, 64 и 128.

inputChannelSize = 67;
hiddenChannelSize = [64,64,128];
nontrainable.nClusters2 = 256;
nontrainable.nSamples2 = 64;
nontrainable.radius2 = 0.2;
trainable.SharedMLP2 = initializeSharedMLP(inputChannelSize,hiddenChannelSize);

Установите параметры третьего модуля абстракции набора. Установите количество кластеров равное 64, количество выборок в каждом кластере равное 32 и радиус каждого кластера равный 0,4. Также установите размер входного канала модели PointNet равным 131, а скрытые размеры каналов равными 128, 128 и 256.

inputChannelSize = 131;
hiddenChannelSize = [128,128,256];
nontrainable.nClusters3 = 64;
nontrainable.nSamples3 = 32;
nontrainable.radius3 = 0.4;
trainable.SharedMLP3 = initializeSharedMLP(inputChannelSize,hiddenChannelSize);

Установите параметры четвертого модуля абстракции набора. Установите количество кластеров равное 16, количество выборок в каждом кластере равное 32 и радиус каждого кластера равный 0,8. Также установите размер входного канала модели PointNet равным 259, а скрытые размеры каналов равными 256, 256 и 512.

inputChannelSize = 259;
hiddenChannelSize = [256,256,512];
nontrainable.nClusters4 = 16;
nontrainable.nSamples4 = 32;
nontrainable.radius4 = 0.8;
trainable.SharedMLP4 = initializeSharedMLP(inputChannelSize,hiddenChannelSize);

Установите параметры первого модуля распространения функций. Установите размер входного канала пятой общей модели MLP равный 768, а размер скрытого канала равный 256 и 256.

inputChannelSize = 768;
hiddenChannelSize = [256,256];
trainable.SharedMLP5 = initializeSharedMLP(inputChannelSize,hiddenChannelSize);

Установите параметры второго модуля распространения функций. Установите размер входного канала общей модели MLP равным 384, а размер скрытого канала равный 256 и 256.

inputChannelSize = 384;
hiddenChannelSize = [256,256];
trainable.SharedMLP6 = initializeSharedMLP(inputChannelSize,hiddenChannelSize);

Установите параметры третьего модуля распространения функций. Установите размер входного канала общей модели MLP равным 320, а размер скрытого канала равный 256 и 128.

inputChannelSize = 320;
hiddenChannelSize = [256,128];
trainable.SharedMLP7 = initializeSharedMLP(inputChannelSize,hiddenChannelSize);

Установите параметры четвертого модуля распространения функций. Установите размер входного канала общей модели MLP равным 128, а размер скрытого канала равным 128 и 128.

inputChannelSize = 128;
hiddenChannelSize = [128,128];
trainable.SharedMLP8 = initializeSharedMLP(inputChannelSize,hiddenChannelSize);

Установите размер входного канала конечной общей модели MLP равным 128, а размер скрытого канала равным 128 и numClasses.

inputChannelSize = 128;
hiddenChannelSize = [128,numClasses];
trainable.SharedMLP9 = initializeSharedMLP(inputChannelSize,hiddenChannelSize);
params.trainable = trainable;
params.nontrainable = nontrainable;

Настройка опций обучения

Обучайте на одну эпоху. Установите начальную скорость обучения равный 2e-4, а L2 коэффициент регуляризации равным 0,01.

numEpochs = 1;
miniBatchSize = 1;
learningRate = 2e-4;
l2Regularization = 0.01;

Инициализируйте опции оптимизации Адама.

gradientDecayFactor = 0.9;
squaredGradientDecayFactor = 0.999;

Инициализируйте скользящее среднее значение градиентов параметра и поэлементных квадратов градиентов, используемых оптимизатором Адама.

trailingAvg = [];
trailingAvgSq = [];

Обучите модель

Обучите модель с помощью CPU или GPU. Для использования GPU требуется Parallel Computing Toolbox™ и графический процессор с поддержкой CUDA ® NVIDIA ®. Для получения дополнительной информации смотрите Поддержку GPU by Release (Parallel Computing Toolbox). Чтобы автоматически обнаружить, доступен ли вам графический процессор, установите executionEnvironment на 'auto'. Если у вас нет графический процессор или вы не хотите использовать его для обучения, установите executionEnvironment на 'cpu'. Чтобы гарантировать использование графический процессор для обучения, установите executionEnvironment на 'gpu'.

Далее создайте minibatchqueue (Deep Learning Toolbox) объект для загрузки данных пакетами miniBatchSize во время обучения.

executionEnvironment = "auto";
if canUseParallelPool
    dispatchInBackground = true;
else
    dispatchInBackground = false;
end

mbq = minibatchqueue(dsTrain,2,...
                     "MiniBatchSize",miniBatchSize,...
                     "OutputEnvironment",executionEnvironment,...
                     "MiniBatchFormat",["SSCB","SSCB"],...
                     "DispatchInBackground",dispatchInBackground);

Обучите модель с помощью пользовательского цикла обучения. Для каждой итерации:

Считайте пакет данных.
Оцените градиенты модели.
Применить L2 регуляризацию веса.
Использование adamupdate для обновления параметров модели.
Обновите график процесса обучения с помощью функции helper configureTrainingProgressPlotter, заданный в конце этого примера.

doTraining = false;
if doTraining
 
    % Initialize plot.
    fig = figure;
    lossPlotter = configureTrainingProgressPlotter(fig);    
    iteration = 0;
    
    % Custom training loop.
    for epoch = 1:numEpochs
         
        % Reset datastore.
        reset(mbq);
        
        while(hasdata(mbq))
            iteration = iteration + 1;
            
            % Read batch of data.
            [ptCld, ptLabels] = next(mbq);
                        
            % Evaluate the model gradients and loss using dlfeval and the modelGradients function.
            [gradients, loss] = dlfeval(@modelGradients, ptCld, ptLabels, params);
                                        
            % L2 regularization.
            gradients = dlupdate(@(g,p) g + l2Regularization*p,gradients,params.trainable);
            
            % Update the network parameters using the Adam optimizer.
            [params.trainable, trailingAvg, trailingAvgSq] = adamupdate(params.trainable, gradients, ...
                                                       trailingAvg, trailingAvgSq, iteration,...
                                                       learningRate,gradientDecayFactor, squaredGradientDecayFactor);
            
            % Update training plot with new points.         
            addpoints(lossPlotter, iteration,double(gather(extractdata(loss))));
            title("Training Epoch " + epoch +" of " + numEpochs);
            drawnow;
        end      
    end
else
    % Download pretrained model parameters.
    pretrainedNetURL = 'https://www.mathworks.com/supportfiles/lidar/data/trainedPointNetplusplusNet.zip';
    
    preTrainedZipFile = fullfile(dataFolder,'trainedPointNetplusplusNet.zip');
    preTrainedMATFile = fullfile(dataFolder,'trainedPointNetplusplusNet.mat');
     if ~exist(preTrainedMATFile,'file')
            if ~exist(preTrainedZipFile,'file')
                disp('Downloading pretrained detector (3.2 MB)...');
                websave(preTrainedZipFile, pretrainedNetURL);
            end
            unzip(preTrainedZipFile,dataFolder);   
     end
    
    % Load pretrained model.
    load(preTrainedMATFile,'params');
    
    % Move model parameters to the GPU if possible and convert to a dlarray.
    params.trainable = prepareForPrediction(params.trainable,@(x)dlarray(toDevice(x,canUseGPU)));   
end

Сегмент Aerial Облако Точек

Ознакомьтесь с полным облаком тестовых точек.

lasReader = lasFileReader(fullfile(testDataFolder,'5175_54395.las'));
[pc,attr] = readPointCloud(lasReader,"Attributes", "Classification");

Вычислим количество непересекающихся сеток на основе gridSize, XLimits, и YLimits облака точек.

numGridsX = round(diff(pc.XLimits)/gridSize(1));
numGridsY = round(diff(pc.YLimits)/gridSize(2));
xRange = linspace(pc.XLimits(1),pc.XLimits(2),numGridsX+1);
yRange = linspace(pc.YLimits(1),pc.YLimits(2),numGridsY+1);

Выполнить итерацию по всем непересекающимся сеткам и предсказать метки.

pcFinal = [];
for r = 1:numGridsX
    for c = 1:numGridsY
        idx = (pc.Location(:,1) >= xRange(r)) & (pc.Location(:,1) < xRange(r+1))...
             & (pc.Location(:,2) >= yRange(c)) & (pc.Location(:,2) < yRange(c+1));
        ptCloudGrid = select(pc,idx);
        
        % Select fixed number of points and labels from the given point cloud.
        ptCloudGrid = selectPoints(ptCloudGrid,numPoints);
        
        % Apply the preprocessing.
        ptCloud = normalizePointCloud(ptCloudGrid);
        loc = ptCloud{1,1};
        X1 = dlarray(loc,'SSCB');
        if (executionEnvironment == "auto" && canUseGPU) || executionEnvironment == "gpu"
            X1 = gpuArray(X1);
        end
        
        % Get the output predictions.
        op = pointNetplusModel(X1,params.trainable,params.nontrainable);
        [~,opLabels] = max(gather(extractdata(op)),[],3);
        color = cmap(opLabels,:);
        ptCloudGrid.Color = uint8(color); 
        
        pcFinal = [pcFinal;ptCloudGrid];
    end
end

ptCloudOut = pccat(pcFinal);
color = 255.*ptCloudOut.Color;
figure;
pcshow(ptCloudOut.Location,color);
title('PointCloud overlaid with detected semantic labels');

Оценка сети

Оцените эффективность сети по тестовым данным. Используйте evaluateSemanticSegmentation функция для вычисления метрики семантической сегментации из результатов тестового набора.

confusionMatrix = {};

% Define the number of samples to evaluate the network. Set numSamplesToTest to
% 1100 to evaluate the model on the entire test data set.
numSamplesToTest = 50;

for i = 1:numSamplesToTest
    ptCldPath = fullfile(pcCropTestPath,sprintf("%03d.pcd",i));
    ptCloud = pcread(ptCldPath);
    ptCloud = normalizePointCloud(ptCloud);
    loc = ptCloud{1,1};
    
    X1 = dlarray(loc,'SSCB');
    if (executionEnvironment == "auto" && canUseGPU) || executionEnvironment == "gpu"
        X1 = gpuArray(X1);
    end
    
    % Get the output predictions.
    op = pointNetplusModel(X1,params.trainable,params.nontrainable);
    [~,predictionLabels] = max(gather(extractdata(op)),[],3);
  
    targetLabelPath = fullfile(labelsCropTestPath,sprintf("%03d.png",i));
    targetLabels = imread(targetLabelPath);
    targetLabels = squeeze(targetLabels);
    
    % Calculate the confusion matrix.
    confMatrix = segmentationConfusionMatrix(double(predictionLabels),...
                 double(targetLabels),"Classes",1:numClasses);
    if size(confMatrix,1) ~= numel(classNames)
        continue
    end
    confusionMatrix{i} = confMatrix;
end

confusionMatrix = confusionMatrix';
metrics = evaluateSemanticSegmentation(confusionMatrix,classNames)

Evaluating semantic segmentation results
----------------------------------------
* Selected metrics: global accuracy, class accuracy, IoU, weighted IoU.
* Processed 50 images.
* Finalizing... Done.
* Data set metrics:

    GlobalAccuracy    MeanAccuracy    MeanIoU    WeightedIoU
    ______________    ____________    _______    ___________

       0.76168          0.53687       0.42486      0.60851

metrics = 
  semanticSegmentationMetrics with properties:

              ConfusionMatrix: [4×4 table]
    NormalizedConfusionMatrix: [4×4 table]
               DataSetMetrics: [1×4 table]
                 ClassMetrics: [4×2 table]
                 ImageMetrics: [50×4 table]

Моделирование градиентов

The modelGradients функция принимает за вход мини-пакет данных X, соответствующий целевой yTarget, и настраиваемые параметры, и возвращает градиенты потери относительно настраиваемых параметров и соответствующих потерь. Чтобы вычислить градиенты, вычислите modelGradients функция, использующая функцию dlfeval в цикле обучения.

function [gradients, loss] = modelGradients(X,yTarget,params)

    % Execute the model function.
    yPred = pointNetplusModel(X,params.trainable,params.nontrainable);
    
    % Compute the loss.
    loss = focalCrossEntropy(yPred,yTarget,'TargetCategories','independent');
    
    % Compute the parameter gradients with respect to the loss. 
    gradients = dlgradient(loss, params.trainable);  
end

Функция модели PointNet++

Функция model принимает как вход параметры модели parameters and входные данные dlX. Сеть выводит предсказания для меток.

function dlY = pointNetplusModel(dlX,trainable,nontrainable)

    % Set abstraction module 1
   [points1, pointFeatures1] = sampleAndGroup(dlX,[],nontrainable.nClusters1,...
                                nontrainable.nSamples1,nontrainable.radius1);
    dlY = sharedMLP(pointFeatures1,trainable.SharedMLP1.Perceptron);
    pointNetFeatures1 = max(dlY,[],2);
    
    % Set abstraction module 2
    [points2, pointFeatures2] = sampleAndGroup(points1,pointNetFeatures1,nontrainable.nClusters2,...
                                nontrainable.nSamples2,nontrainable.radius2);
    dlY = sharedMLP(pointFeatures2,trainable.SharedMLP2.Perceptron);
    pointNetFeatures2 = max(dlY,[],2);
    
    % Set abstraction module 3
    [points3, pointFeatures3] = sampleAndGroup(points2,pointNetFeatures2,nontrainable.nClusters3,...
                                nontrainable.nSamples3,nontrainable.radius3);  
    dlY = sharedMLP(pointFeatures3,trainable.SharedMLP3.Perceptron);
    pointNetFeatures3 = max(dlY,[],2);
    
    % Set abstraction module 4
   [points4, pointFeatures4] = sampleAndGroup(points3,pointNetFeatures3,nontrainable.nClusters4,...
                                nontrainable.nSamples4,nontrainable.radius4); 
    dlY = sharedMLP(pointFeatures4,trainable.SharedMLP4.Perceptron);
    pointNetFeatures4 = max(dlY,[],2);
    
    % Feature propagation module 1
    pointsFP1 = featurePropagation(points3,points4,pointNetFeatures3,pointNetFeatures4);
    pointNetFP1 = sharedMLP(pointsFP1,trainable.SharedMLP5.Perceptron);
    
    % Feature propagation module 2
    pointsFP2 = featurePropagation(points2,points3,pointNetFeatures2,pointNetFP1);
    pointNetFP2 = sharedMLP(pointsFP2,trainable.SharedMLP6.Perceptron);
   
    % Feature propagation module 3
    pointsFP3 = featurePropagation(points1,points2,pointNetFeatures1,pointNetFP2);
    pointNetFP3 = sharedMLP(pointsFP3,trainable.SharedMLP7.Perceptron);
    
    % Feature propagation module 4
    pointsFP4 = featurePropagation(dlX,points1,[],pointNetFP3);
    dlY = sharedMLP(pointsFP4,trainable.SharedMLP8.Perceptron);
   
    % Shared MLP
    dlY = sharedMLP(dlY,trainable.SharedMLP9.Perceptron);
    dlY = softmax(dlY);
    dlY = dlarray(dlY,'SSCB');
end

Функции инициализации параметра модели

Инициализируйте общую многослойную функцию перцептрона

The initializeSharedMLP функция принимает за вход размер канала и скрытый размер канала, и возвращает инициализированные параметры в структуре. Параметры инициализируются с помощью инициализации веса He.

function params = initializeSharedMLP(inputChannelSize,hiddenChannelSize)
weights = initializeWeightsHe([1 1 inputChannelSize hiddenChannelSize(1)]);
bias = zeros(hiddenChannelSize(1),1,"single");
p.Conv.Weights = dlarray(weights);
p.Conv.Bias = dlarray(bias);

p.InstanceNorm.Offset = dlarray(zeros(hiddenChannelSize(1),1,"single"));
p.InstanceNorm.Scale = dlarray(ones(hiddenChannelSize(1),1,"single"));

params.Perceptron(1) = p;

for k = 2:numel(hiddenChannelSize)
    weights = initializeWeightsHe([1 1 hiddenChannelSize(k-1) hiddenChannelSize(k)]);
    bias = zeros(hiddenChannelSize(k),1,"single");
    p.Conv.Weights = dlarray(weights);
    p.Conv.Bias = dlarray(bias);
    
    p.InstanceNorm.Offset = dlarray(zeros(hiddenChannelSize(k),1,"single"));
    p.InstanceNorm.Scale = dlarray(ones(hiddenChannelSize(k),1,"single"));
   
    params.Perceptron(k) = p;
end
end

function x = initializeWeightsHe(sz)
fanIn = prod(sz(1:3));
stddev = sqrt(2/fanIn);
x = stddev .* randn(sz);
end

function x = initializeWeightsGaussian(sz)
x = randn(sz,"single") .* 0.01;
end

Выборка и группировка

function [newClusters,newpointFeatures] = sampleAndGroup(points,pointFeatures,nClusters,nSamples,radius)
% The sampleAndGroup layer first samples the point cloud to a given number of
% clusters and then constructs local region sets by finding neighboring
% points around the centroids using the queryBallPoint function.

    points = extractdata(squeeze(points));
    if ~isempty(pointFeatures)
        pointFeatures = extractdata(squeeze(pointFeatures));
    end
        
    % Find the centroids for nClusters - nClusters*3.
    centroids = farthestPointSampling(points,nClusters);
    newClusters = points(centroids,:);
    
    % Find the neareset nSamples for nClusters - nClusters*nSamples*3.
    idx = queryBallPoint(points,newClusters,nClusters,nSamples,radius);
    newPoints = reshape(points(idx,:),[nClusters,nSamples,3]);
    
    % Normalize the points in relation to the cluster center.
    newpointFeatures = newPoints - reshape(newClusters,nClusters,1,3);
    
    if ~isempty(pointFeatures)
        groupFeatures = reshape(pointFeatures(idx,:),...
                       [nClusters,nSamples,size(pointFeatures,2)]);
        newpointFeatures = cat(3,newPoints,groupFeatures);
    end
    
    newpointFeatures = dlarray(newpointFeatures,'SSC');
    newClusters = dlarray(newClusters,'SSC');
end

Отбор проб в самой дальней точке

function centroids = farthestPointSampling(pointLocations,numPoints)
% The farthestPointSampling function selects a set of points from input
% points, which defines the centroids of local regions.
% pointLocations - PointCloud locations N-by-3.
% numPoints - Number of clusters to find.
% centroids - centroids of each farthest cluster.
    
    % Initialize initial indices as zeros.
    centroids = zeros(numPoints,1);
    
    % Distance from centroid to each point.
    distance = ones(size(pointLocations,1),1) .* 1e10; 
    
    % Random Initialization of the first point.
    farthest = randi([1,size(pointLocations,1)],1);
    
    for i = 1:numPoints
        centroids(i) = farthest;
        centroid = pointLocations(farthest,:);
        dist = sum(((pointLocations - centroid).^2),2);
        mask = dist < distance;
        distance(mask) = dist(mask);
        [~,farthest] = max(distance,[],1);
    end
end

Запрос шаровой точки

function groupIdx = queryBallPoint(XYZ,newXYZ,nClusters,nSamples,radius)
% Given the cluster center, the queryBallPoint finds all points that are
% within a radius to the query point.

    N = size(XYZ,1);
    groupIdx = reshape(1:N,[1,N]);
    groupIdx = repmat(groupIdx,[nClusters,1]);
    
    % Find the distance between the centroids and given points.
    sqDist = squareDistance(newXYZ,XYZ);    
    
    % Find the points that are inside the given radius.
    groupIdx(sqDist > (radius)^2) = N;
    groupIdx = sort(groupIdx,2,"ascend");
    
    % Find the closest nSamples points within the given radius.
    groupIdx = groupIdx(:,1:nSamples);
    groupFirst = repmat(groupIdx(:,1),1,nSamples);
    mask = (groupIdx == N);
    groupIdx(mask) = groupFirst(mask);
end

Распространение функций

function newPoints = featurePropagation(points1,points2,pointNetFeatures1,pointNetFeatures2)
    % Use the inverse distance weighted average based on the k nearest neighbors to
    % interpolate features.
    
    points1 = extractdata(squeeze(points1));
    points2 = extractdata(squeeze(points2));
    if ~isempty(pointNetFeatures1)
    pointNetFeatures1 = extractdata(squeeze(pointNetFeatures1));
    end
    pointNetFeatures2 = extractdata(squeeze(pointNetFeatures2));
    
    % Find the K nearest neighbors for each point.
    dists = squareDistance(points1,points2);
    [dists,idx] = sort(dists,2,"ascend");
    dists = dists(:,1:3);
    idx = idx(:,1:3);
    
    % Calculate the weights for interpolation.
    dist_recip = 1./(dists+1e-8);
    normFactor = sum(dist_recip,2);
    weights = dist_recip./normFactor;
    
    % Perform weighted interpolation.
    interpolatedPoints = pointNetFeatures2(idx,:);
    interpolatedPoints = reshape(interpolatedPoints,[size(idx,1),size(idx,2),size(pointNetFeatures2,2)]);
    
    interpolatedPoints = interpolatedPoints .* weights;
    interpolatedPoints = squeeze(sum(interpolatedPoints,2));
    
    if ~isempty(pointNetFeatures1)
        % Calculate the new points.
        newPoints = cat(2,pointNetFeatures1,interpolatedPoints);
    else
        newPoints = interpolatedPoints;
    end
    
     newPoints = dlarray(newPoints,'SCS');
end

% Find the squared distance.
function dist = squareDistance(src,dst)
    dist = -2 * (src*permute(dst,[2,1,3]));
    tmp1 = sum(src.^2,2);
    tmp1 = reshape(tmp1,[size(src,1),1]);
    tmp2 = sum(dst.^2,2);
    tmp2 = reshape(tmp2,[1,size(dst,1)]);
    dist = dist + tmp1 + tmp2;
end

Общая многослойная функция перцептрона

function dlY= sharedMLP(dlX,parameters)
% The shared multilayer perceptron function processes the input dlX using a
% series of perceptron operations and returns the result of the last
% perceptron.
    dlY = dlX;
    for k = 1:numel(parameters) 
        dlY = perceptron(dlY,parameters(k));
    end
end

Функция Перцептрона

function dlY = perceptron(dlX,parameters)
% The perceptron function processes the input dlX using a convolution, a
% instance normalization, and a ReLU operation and returns the output of the
% ReLU operation.

    % Convolution.
    W = parameters.Conv.Weights;
    B = parameters.Conv.Bias;
    dlY = dlconv(dlX,W,B);
    
    % Instance normalization.
    offset = parameters.InstanceNorm.Offset;
    scale = parameters.InstanceNorm.Scale;
    dlY = instancenorm(dlY,offset,scale);
    
    % ReLU.
    dlY = relu(dlY);
end

Вспомогательные функции

The normalizePointCloud функция извлекает данные точек X, Y, Z из входных данных и нормализует данные между 0 и 1. Функция возвращает нормированные данные X, Y, Z.

function data = normalizePointCloud(data)
    if ~iscell(data)
        data = {data};
    end
    
    numObservations = size(data,1);
    for i = 1:numObservations
        % Scale points between 0 and 1.
        xlim = data{i,1}.XLimits;
        ylim = data{i,1}.YLimits;
        zlim = data{i,1}.ZLimits;
        
        xyzMin = [xlim(1) ylim(1) zlim(1)];
        xyzDiff = [diff(xlim) diff(ylim) diff(zlim)];
        
        tmp = (data{i,1}.Location - xyzMin) ./ xyzDiff;
        % Convert the data to SCSB format.
        data{i,1} = permute(tmp,[1,3,2]);
    end
end

function data = onehotEncodeLabels(data,classNames)
    numObservations = size(data,1);
    for i = 1:numObservations
       labels = data{i,1}';
       encodedLabels = onehotencode(labels,2,'ClassNames',classNames);
       data{i,1} = permute(encodedLabels,[1,3,2]);
    end
end

`prepareForPrediction` Функция

The prepareForPrediction функция используется для применения пользовательской функции к данным вложенной структуры. Используйте эту функцию, чтобы переместить параметр модели и данные о состоянии в графический процессор.

function p = prepareForPrediction(p,fcn)

for i = 1:numel(p)
    p(i) = structfun(@(x)invoke(fcn,x),p(i),'UniformOutput',0);
end

    function data = invoke(fcn,data)
        if isstruct(data)
            data = prepareForPrediction(data,fcn);
        else
            data = fcn(data);
        end
    end
end

% Move data to the GPU.
function x = toDevice(x,useGPU)
if useGPU
    x = gpuArray(x);
end
end

`selectPoints` Функция

The selectPoints функция производит выборку нужного числа точек. Когда облако точек содержит больше, чем желаемое число точек, функция использует pcdownsample для случайного выбора точек. В противном случае функция наследует данные для создания требуемого числа точек.

function ptCloudOut = selectPoints(ptCloud,numPoints) 
    if ptCloud.Count > numPoints
        ind = 1:ptCloud.Count;
    else    
        replicationFactor = ceil(numPoints/ptCloud.Count);
        ind = repmat(1:ptCloud.Count,1,replicationFactor);
    end 
    
     ptCloudOut = select(ptCloud,ind(1:numPoints));
end

The labelColorbar функция добавляет шкалу палитры к текущей оси. Шкала палитры форматирована, чтобы отображать имена классов с цветом.

function labelColorbar(ax, cmap, classNames)
    colormap(ax, cmap);
    
    % Add colorbar to current figure.
    c = colorbar(ax);
    c.Color = 'w';
        
    % Center tick labels and use class names for tick marks.
    numClasses = size(classNames, 1);
    c.Ticks = 1:1:numClasses;
    c.TickLabels = classNames; 

    % Remove tick mark.
    c.TickLength = 0;
end

function cmap = labelColorMap()
% Colormap for the original classes.
cmap = [[10,10,10];
        [0,0,255];
        [0,255,0];
        [255,192,203];
        [255,255,0];
        [255,0,255];
        [255,165,0];
        [139,0,150];
        [255,0,0]];
cmap = cmap./255;
end

function lossPlotter = configureTrainingProgressPlotter(f)
% The configureTrainingProgressPlotter function configures training
% progress plots for various losses.
    figure(f);
    clf
    ylabel('Total Loss');
    xlabel('Iteration');
    lossPlotter = animatedline;
end

Ссылки

[1] Варней, Нина, Виджаян К. Асари и Куинн Грелинг. DALES: Крупномасштабный набор данных воздушного ЛИДАРА для семантической сегментации. ArXiv:2004.11985 [Cs, Stat], 14 апреля 2020 года. https://arxiv.org/abs/2004.11985.

[2] Qi, Charles R., Li Yi, Hao Su, and Leonidas J. Guibas. «PointNet++: Глубокое иерархическое обучение функций на наборах точек в метрическом пространстве». ArXiv:1706.02413 [Cs], 7 июня 2017 года. https://arxiv.org/abs/1706.02413.

Документация

Семантическая сегментация Aerial Lidar с использованием глубокого обучения PointNet++

Загрузка данных DALES

Предварительная обработка данных

Создайте объекты Datastore для обучения

Задайте модель PointNet++

Задайте параметры модели PointNet++

Настройка опций обучения

Обучите модель

Сегмент Aerial Облако Точек

Оценка сети

Моделирование градиентов

Функция модели PointNet++

Функции инициализации параметра модели

Инициализируйте общую многослойную функцию перцептрона

Выборка и группировка

Отбор проб в самой дальней точке

Запрос шаровой точки

Распространение функций

Общая многослойная функция перцептрона

Функция Перцептрона

Вспомогательные функции

`prepareForPrediction` Функция

`selectPoints` Функция

Ссылки

Документация по Lidar Toolbox

Поддержка

Документация

Семантическая сегментация Aerial Lidar с использованием глубокого обучения PointNet++

Загрузка данных DALES

Предварительная обработка данных

Создайте объекты Datastore для обучения

Задайте модель PointNet++

Задайте параметры модели PointNet++

Настройка опций обучения

Обучите модель

Сегмент Aerial Облако Точек

Оценка сети

Моделирование градиентов

Функция модели PointNet++

Функции инициализации параметра модели

Инициализируйте общую многослойную функцию перцептрона

Выборка и группировка

Отбор проб в самой дальней точке

Запрос шаровой точки

Распространение функций

Общая многослойная функция перцептрона

Функция Перцептрона

Вспомогательные функции

prepareForPrediction Функция

selectPoints Функция

Ссылки

Документация по Lidar Toolbox

Поддержка

`prepareForPrediction` Функция

`selectPoints` Функция