Main Content

Aerial Lidar Semantic Segmentation Using PointNet++ Deep Learning

This example shows how to train a PointNet++ deep learning network to perform semantic segmentation on aerial lidar data.

Lidar data acquired from airborne laser scanning systems is increasingly used in applications such as topographic mapping, building and city modeling, biomass measurement, and disaster management. Extracting meaningful information from this data requires semantic segmentation, a process where each pixel in the point cloud is assigned a class label.

In this example, you train a PointNet++ network to perform semantic segmentation by using the Dayton annotated lidar earth scan (DALES) data set [1], which contains scenes of dense, labeled lidar aerial data from urban, suburban, rural, and commercial settings. Out of 40 scenes, 29 scenes are used for training and the remaining 11 sequences are used for testing. The data set provides semantic segmentation labels for 8 such as buildings, cars, trucks, poles, power lines, fences, ground, and vegetation.

Load DALES Data

The DALES data set contains 40 scenes of aerial lidar data. The data set is already divided into 29 scenes for training and 11 scenes for testing. Each pixel in the data has a class label. Follow the instructions on the DALES website to download the data set to the folder specified by the dataFolder variable. Create folders to store training and test data.

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

Preview a point cloud from the training data.

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

Preprocess Data

Each point cloud in the DALES data set covers an area of 500-by-500 meters, which is much larger than the typical area covered by terrestrial rotating lidar point clouds. For efficient memory processing, divide the point cloud into small, non-overlapping grids. Additionally, focus on the three dominant classes—the ground, buildings, and vegetation—to avoid problems due to class imbalance.

Use the helperCropPointCloudsAndMergeLabels function, attached to this example as a supporting file, to:

  • Crop the point clouds into non-overlapping grids of size 50-by-50 meters.

  • Merge all classes other than the dominant classes (the ground, buildings, and vegetation) into a single background class so that the data contains only four classes.

  • Save the cropped grids and semantic labels as PCD and PNG files, respectively.

Define the grid dimensions and set a fixed number of points per grid to enable faster training.

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

If you use this workflow for your own data, set writeFiles to "false" if the training data is divided into grids. Please note that the training data must be in a format supported by the pcread function.

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

Processing can take some time. The code suspends MATLAB® execution until processing is complete.

Create Datastore Objects for Training

Create a fileDatastore object to load PCD files using the pcread function.

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

Use a pixelLabelDatastore object to store pixel-wise labels from the pixel label images. The object maps each pixel label to a class name. In this example, the ground, vegetation, and buildings are the only classes of interest; all other pixels are the background. Specify these classes and assign a unique label ID to each class.

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

Load and display the point clouds.

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

Use the normalizePointCloud function, defined at the end of the example, to normalize the point cloud between 0 and 1.

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

Use the onehotEncodeLabels function, defined at the end of the example, to convert the labels to one-hot-encoded representation for training the segmentation network using a custom training loop.

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

Use the combine function to combine the point clouds and pixel labels into a single datastore for training.

dsTrain = combine(ldsTransformed,pxdsTransformed);

Define PointNet++ Model

The PointNet++ [2] segmentation model consists of two main components:

  • Set abstraction modules

  • Feature propagation modules

The series of set abstraction modules progressively subsamples points of interest by hierarchically grouping points, and uses a custom PointNet architecture to encode points into feature vectors. Because semantic segmentation tasks require point features for all original points, this example uses a series of feature propagation modules to hierarchically interpolate features to original points using an inverse- distance-based interpolation scheme.

The set abstraction module is implemented using the sampleAndGroup and sharedMLP supporting functions and the feature propagation module is implemented using the featurePropagation and sharedMLP supporting functions. The sharedMLP module is the only learnable module that consists of convolution and instance normalization layers. Use the supporting function pointNetplusModel to define the entire PointNet++ segmentation model. The supporting functions are listed at the end of this example.

Define PointNet++ Model Parameters

The sampleAndGroup function of the set abstraction module is parameterized by the number of clusters, number of samples per cluster, and radius of each cluster, whereas the sharedMLP function implements a custom PointNet architecture that is parameterized by using the number of input channels and the hidden channel sizes.

This example uses hyperparameter values tuned on the DALES data set. If you want to apply PointNet++ to a different data set, you must perform additional hyperparameter tuning.

Set the parameters of the first set abstraction module. Set the number of clusters to 1024, the number of samples in each cluster to 128, and the radius of each cluster to 0.1. Also set the PointNet model input channel size to three and the hidden channel sizes to 32, 32, and 64.

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

Set the parameters of the second set abstraction module. Set the number of clusters to 256, the number of samples in each cluster to 64, and the radius of each cluster to 0.2. Also set the PointNet model input channel size to 67 and the hidden channel sizes to 64, 64, and 128.

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

Set the parameters of the third set abstraction module. Set the number of clusters to 64, the number of samples in each cluster to 32, and the radius of each cluster to 0.4. Also set the PointNet model input channel size to 131 and the hidden channel sizes to 128, 128, and 256.

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

Set the parameters of the fourth set abstraction module. Set the number of clusters to 16, the number of samples in each cluster to 32, and the radius of each cluster to 0.8. Also set the PointNet model input channel size to 259 and the hidden channel sizes to 256, 256, and 512.

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

Set the parameters of the first feature propagation module. Set the fifth shared MLP model input channel size to 768 and the hidden channel size to 256 and 256.

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

Set the parameters of the second feature propagation module. Set the shared MLP model input channel size to 384 and the hidden channel size to 256 and 256.

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

Set the parameters of the third feature propagation module. Set the shared MLP model input channel size to 320 and the hidden channel size to 256 and 128.

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

Set the parameters of the fourth feature propagation module. Set the shared MLP model input channel size to 128 and the hidden channel size to 128 and 128.

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

Set the final shared MLP model input channel size to 128 and the hidden channel size to 128 and numClasses.

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

Specify Training Options

Train for one epoch. Set the initial learning rate to 2e-4 and the L2 regularization factor to 0.01.

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

Initialize the options for Adam optimization.

gradientDecayFactor = 0.9;
squaredGradientDecayFactor = 0.999;

Initialize the moving average of the parameter gradients and the element-wise squares of the gradients used by the Adam optimizer.

trailingAvg = [];
trailingAvgSq = [];

Train Model

Train the model using a CPU or GPU. Using a GPU requires Parallel Computing Toolbox™ and a CUDA® enabled NVIDIA® GPU. For more information, see GPU Support by Release (Parallel Computing Toolbox). To automatically detect if you have a GPU available, set executionEnvironment to 'auto'. If you do not have a GPU, or do not want to use one for training, set executionEnvironment to 'cpu'. To ensure the use of a GPU for training, set executionEnvironment to 'gpu'.

Next, create a minibatchqueue (Deep Learning Toolbox) object to load the data in batches of miniBatchSize during training.

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

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

Train the model using a custom training loop. For each iteration:

  • Read a batch of data.

  • Evaluate the model gradients.

  • Apply L2 weight regularization.

  • Use adamupdate to update the model parameters.

  • Update the training progress plotter using the helper function configureTrainingProgressPlotter, defined at the end of this example.

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

Segment Aerial Point Cloud

Read the full test point cloud.

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

Calculate the number of non-overlapping grid based on gridSize, XLimits, and YLimits of the point cloud.

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

Iterate over all the non-overlapping grids and predict the labels.

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

Evaluate Network

Evaluate the network performance on the test data. Use the evaluateSemanticSegmentation function to compute the semantic segmentation metrics from the test set results.

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]

Model Gradients

The modelGradients function takes as input a mini-batch of data X, the corresponding target yTarget, and the learnable parameters, and returns the gradients of the loss with respect to the learnable parameters and the corresponding loss. To compute the gradients, evaluate the modelGradients function using the dlfeval function in the training loop.

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 Function

The model function takes as input the model parameters parameters and the input data dlX. The network outputs the predictions for the labels.

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

Model Parameter Initialization Functions

Initialize Shared Multilayer Perceptron Function

The initializeSharedMLP function takes as input the channel size and the hidden channel size, and returns the initialized parameters in a structure. The parameters are initialized using He weight initialization.

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

Sampling and Grouping

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

Farthest Point Sampling

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

Query Ball Point

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

Feature Propagation

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

Shared Multilayer Perceptron Function

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

Perceptron Function

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

Supporting Functions

The normalizePointCloud function extracts the X, Y, Z point data from the input data and normalizes the data between 0 and 1. The function returns the normalized X, Y, Z data.

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 Function

The prepareForPrediction function is used to apply a user-defined function to nested structure data. Use this function to move model parameter and state data to the GPU.

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 Function

The selectPoints function samples the desired number of points. When the point cloud contains more than the desired number of points, the function uses pcdownsample to randomly select points. Otherwise, the function replicates data to produce the desired number of points.

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 adds a colorbar to the current axis. The colorbar is formatted to display the class names with the color.

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

References

[1] Varney, Nina, Vijayan K. Asari, and Quinn Graehling. "DALES: A Large-Scale Aerial LiDAR Data Set for Semantic Segmentation." ArXiv:2004.11985 [Cs, Stat], April 14, 2020. https://arxiv.org/abs/2004.11985.

[2] Qi, Charles R., Li Yi, Hao Su, and Leonidas J. Guibas. "PointNet++: Deep Hierarchical Feature Learning on Point Sets in a Metric Space." ArXiv:1706.02413 [Cs], June 7, 2017. https://arxiv.org/abs/1706.02413.