Determine Safe Landing Area for Aerial Vehicles

Open Live Script

This example demonstrates how to determine a safe landing area for aerial vehicles, such as helicopters and UAVs, by using aerial lidar data.

Aerial vehiclesare increasingly used for applications like cargo delivery, casualty evacuation, and search and rescue operations. For these applications, the aerial vehicle must land safely in the vicinity of the destination. Because landing an aerial vehicle based on limited visual information is a challenging task for a pilot, aerial vehicles benefit from using lidar data to perform, reducing the pressure on the pilot and the risk of an accident.

Use these steps to detect safe landing zones from aerial lidar data:

Load the point cloud into the workspace.
Divide the point cloud into overlapping blocks.
Identify the non-flat points, and label these points as dangerous.
Identify and label the points around the dangerous points as unsuitable.
Identify and label the points around the unsuitable points as risky. Also, label those points that do not satisfy the parameter thresholds as risky.
Identify the safe landing zones.

You must avoid the dangerous and unsuitable points for landing, as these carry a very high chance of causing an accident. The probability of accidents is lower at the risky points, but the terrain points labeled as suitable are the safest for landing.

Load and Visualize Data

This example uses the point cloud data from the LAZ file aerialLidarData.laz, obtained from the Open Topography Dataset [1]. Load the point cloud data into the workspace using the readPointCloud function of the lasFileReader object, and visualize the point cloud.

% Set random seed to generate reproducible results
rng("default")

% Specify the LAZ file path
lasfile = fullfile(toolboxdir("lidar"),"lidardata","las","aerialLidarData.laz");

% Create a lasFileReader object
lasReader = lasFileReader(lasfile);

% Read the point cloud
ptCloud = readPointCloud(lasReader);

% Visualize the point cloud
figure
pcshow(ptCloud)
axis on
title("Input Aerial Lidar Data")

Figure contains an axes object. The axes object with title Input Aerial Lidar Data contains an object of type scatter.

Divide Point Cloud into Blocks

Divide the input point cloud into overlapping blocks. Each block consists of an inner block and an outer block. First, specify the size of the inner block.

% Specify the inner block size (in meters)
innerBlockSize = [20 20];

The outer block size is the sum of the inner block size and landing site radius. Specify the landing site radius.

% Specify the radius of the landing site to be evaluated (in meters)
radius = 5;

The landing site radius is the sum of the airframe radius of the aerial vehicle and the landing error, as shown in this figure.

While labeling, you label only the points in the inner block. You determine the label of each point in the inner block by evaluating the parameters of that point and its nearest neighbors within the landing site radius.

For example, to assign a label to the red point in this figure, you must evaluate the properties of the red point and its neighbors within the radius, represented by the green points.

The block processing starts from the eft and proceeds to the right, then repeats from bottom to top until all the blocks have been covered. Define the overlap between two adjacent blocks such that the inner blocks of the adjacent blocks lay side-by-side without any overlap with each other, as shown in this figure.

Use the helperOverlappingBlocks helper function, attached to this example as a supporting file, to compute the parameters required for overlapping-block-based processing. These parameters contain the indices of the outer blocks, inner blocks, and boundary points. The function also outputs the mapping for each block between the inner block, outer block, and labels.

[outerBlockIndices,innerBlockIndices, ...
    innerBlockToOuterBlock,innerBlockIndFromLabels, ...
    innerBlockIndFromOuterBlock,boundaryIndices] = helperOverlappingBlocks(ptCloud, ...
    radius,innerBlockSize);

Define the class labels for classifying each point.

classNames = [
    "Unlabeled"
    "Dangerous"     % Label 1
    "Unsuitable"    % Label 2
    "Risky"         % Label 3
    "Suitable"      % Label 4
    ];

You use the overlapping block processing parameters when labeling the unsuitable and risky points. Use overlapping block processing to improve the run-time performance of the labeling algorithm on larger point clouds.

Classify Dangerous Points

A good landing zone must be a flat terrain surface with minimal obstacles, as landing an aerial vehicle on non-flat points is dangerous, and can result in an accident. You must further evaluate the flat terrain points to analyze their safety.

Flat terrain points generally consist of the ground points, while non-flat points generally consist of trees, buildings, and other non-ground points.

Segment the input point cloud into ground and non-ground points, using the segmentGroundSMRF function. Label the non-ground points as dangerous.

% Store the label of each point in a labels array
labels = nan(ptCloud.Count,1);

% Identify the ground points
flatIdx = segmentGroundSMRF(ptCloud);

% Mark the non-ground points as dangerous
labels(~flatIdx) = 1;

% Use the helperViewOutput helper function, defined at the end of the
% example, to visualize the output
helperViewOutput(ptCloud,labels,classNames)
title("Dangerous Points")

Figure contains an axes object. The axes object with title Dangerous Points contains an object of type scatter.

Classify Unsuitable Points

Label a point as unsuitable in any of these cases:

The point is along the boundary of the point cloud.
There are dangerous points in the neighborhood of the point.
The point does not have enough neighboring points, and further evaluation is not possible.

Classify Unsuitable Boundary Points

Because you cannot assess the entire neighborhood of the points along the boundary of the point cloud, label any unlabeled boundary points as unsuitable.

% Identify the unlabeled boundary points
unlabeledBoundaryIndices = isnan(labels(boundaryIndices));

% Label these points as unsuitable
labels(boundaryIndices(unlabeledBoundaryIndices)) = 2;

% Use the helperViewOutput helper function, defined at the end of the
% example, to visualize the output
helperViewOutput(ptCloud,labels,classNames)
title("Unsuitable Points Along the Boundary")

Figure contains an axes object. The axes object with title Unsuitable Points Along the Boundary contains an object of type scatter.

Classify Unsuitable Points Around Dangerous Points

Landing an aerial vehicle lands near dangerous points can result in an accident during landing. Dangerous points consist of powerlines, trees, buildings, and other non-ground objects. Due to their height, these points might not be detected as a neighbors within the radius of the flat points. Use the helperUnsuitablePoints helper function to label points in the vicinity of the dangerous points as unsuitable.

% Define the grid size (in meters)
gridSize = 4;

Note: When you increase the gridSize value, the helperUnsuitablePoints function detects more points in the vicinity of each dangerous point, thus labeling more points as unsuitable.

% Identify the indices of the points near the dangerous points
unsuitableIdx = helperUnsuitablePoints(ptCloud,labels,gridSize);

% Label these points as unsuitable
labels(unsuitableIdx) = 2;

% Use the helperViewOutput helper function, defined at the end of the
% example, to visualize the output
helperViewOutput(ptCloud,labels,classNames)
title("Unsuitable Points Near the Dangerous Points")

Figure contains an axes object. The axes object with title Unsuitable Points Near the Dangerous Points contains an object of type scatter.

Extract Indices of Nearest Neighbors for Unlabeled Points

Further classify unsuitable points by using overlapping-block-based processing, and compute the nearest neighbors for each unlabeled point.

If the number of neighbors is less than the minPoints value, or the neighborhood of the point contains dangerous points, label the point as unsuitable. Otherwise, store the indices of the neighboring points in the totalNeighborInd cell array.

% Define the minimum number of neighbors for each point
minPoints = 10;

% Store the indices of the nearest neighboring points around each point
% within the specified radius
totalNeighborInd = {};

% Perform block-by-block processing
for curBlockIdx = 1:numel(outerBlockIndices)
    % Extract the inner block point cloud from the input point cloud
    innerBlockPtCloud = select(ptCloud,innerBlockIndices{curBlockIdx});

    % Extract the outer block point cloud from the input point cloud
    outerBlockPtCloud = select(ptCloud,outerBlockIndices{curBlockIdx});

    % Extract the labels of the outer block
    curOuterBlockLabels = labels(outerBlockIndices{curBlockIdx});

    % Create a mapping from the inner block labels to the outer block labels
    curInnerBlockToOuterBlock = innerBlockToOuterBlock{curBlockIdx};

Use the helperNeighborIndices helper function, defined at the end of the example, to compute the nearest neighbor indices within the radius of each point. The function also labels a point as unsuitable if the number of nearest neighbors within the specified radius of that point is less than the minimum number of points.

    [labeledOuterBlock,neighborInd] = helperNeighborIndices( ...
        innerBlockPtCloud,outerBlockPtCloud,curOuterBlockLabels, ...
        curInnerBlockToOuterBlock,minPoints,radius);

    % Store the neighbor indices of the points belonging to the inner block
    totalNeighborInd = [totalNeighborInd; neighborInd];

Label only the inner block point cloud. Compute the indices of the inner block with respect to the labels and the outer block. Update the labels array with the labels computed for the inner block.

    labels(innerBlockIndFromLabels{curBlockIdx}) = labeledOuterBlock( ...
        innerBlockIndFromOuterBlock{curBlockIdx});
end

% Use the helperViewOutput helper function, defined at the end of the
% example, to visualize the output
helperViewOutput(ptCloud,labels,classNames)
title("Unsuitable Points Detected")

Figure contains an axes object. The axes object with title Unsuitable Points Detected contains an object of type scatter.

Classify Risky Points

Within the unlabeled points, label a point as risky, if it has unsuitable points in its neighborhood.

Also, evaluate these attributes of each unlabeled point and its neighboring points.

Vertical Variance — Variance of the points along the z-axis. A higher variance value indicates a greater height spread among the points, which can make them unsuitable for landing.
Relief — Height difference between the lowest and highest points in the neighborhood. A smaller relief value correlates with a flatter surface and fewer obstacles in the landing zone.
Slope — Inclination angle of the fitting plane. A smaller slope value is more suitable for landing, as indicates a more stable surface for the vehicle.
Residual — Roughness of the landing zone. A smaller roughness value indicates the presence of fewer obstacles and a smoother landing zone.

Note: Because a landing surface must be flat, lower values of the vertical variance and relief attributes ensure that you can fit a plane through the points. Based on the plane fitting, the algorithm computes the slope and residual attributes.

Specify the threshold values for these attributes. Any point, with an attribute value is greater than the respective threshold is a risky point.

% Define the vertical variance threshold (in meters)
verticalVarianceThreshold = 0.5;

% Define the slope threshold (in degrees)
slopeThreshold = 6;

% Define the residual threshold (in meters)
residualThreshold = 10;

% Define the relief threshold (in meters)
reliefThreshold = 3;

innerBlockIndicesCount = 1;
for curBlockIdx = 1:numel(outerBlockIndices)
    % Extract the inner block point cloud from the input point cloud
    innerBlockPtCloud = select(ptCloud,innerBlockIndices{curBlockIdx});

    % Extract the outer block point cloud from the input point cloud
    outerBlockPtCloud = select(ptCloud,outerBlockIndices{curBlockIdx});

    % Extract the labels of the outer block
    curOuterBlockLabels = labels(outerBlockIndices{curBlockIdx});

    % Map the inner block labels to the outer block labels
    curInnerBlockToOuterBlock = innerBlockToOuterBlock{curBlockIdx};

Use the helperLabelRiskyPts helper function, defined at the end of the example, to label the points in the vicinity of the unsuitable points as risky. Also, label the points whose slope, vertical variance, relief, or residual parameter is greater than the corresponding threshold as risky.

    [labeledOuterBlock,updatedInnerBlockIndicesCount, ...
        updatedTotalNeighborInd] = helperLabelRiskyPts( ...
        innerBlockIndicesCount,totalNeighborInd,innerBlockPtCloud, ...
        outerBlockPtCloud,curOuterBlockLabels, ...
        curInnerBlockToOuterBlock, ...
        verticalVarianceThreshold,slopeThreshold,residualThreshold,...
        reliefThreshold);

    totalNeighborInd = updatedTotalNeighborInd;
    innerBlockIndicesCount = updatedInnerBlockIndicesCount;

Label only the inner block point cloud. Compute the indices of the inner block with respect to the labels and the outer block. Update the labels array with the labels computed for the inner block.

    labels(innerBlockIndFromLabels{curBlockIdx}) = labeledOuterBlock( ...
        innerBlockIndFromOuterBlock{curBlockIdx});
end

% Use the helperViewOutput helper function, defined at the end of the
% example, to visualize the output
helperViewOutput(ptCloud,labels,classNames)
title("Risky Points Detected")

Figure contains an axes object. The axes object with title Risky Points Detected contains an object of type scatter.

Classify Safe Landing Zones

The remaining unlabeled points, which satisfy the parameter threshold values, are the safe landing points. There are no dangerous or unsuitable points around these points. Label these points as suitable.

% Label the remaining unlabeled points as suitable
labels(isnan(labels)) = 4;

% Use the helperViewOutput helper function, defined at the end of the
% example, to visualize the output
helperViewOutput(ptCloud,labels,classNames)
title("Landing Zones Detected")

Figure contains an axes object. The axes object with title Landing Zones Detected contains an object of type scatter.

Helper Functions

The helperViewOutput helper function visualizes the point cloud along with the class labels.

function helperViewOutput(ptCloud,labels,classNames)
labels = labels + 1;
labels(isnan(labels)) = 1;

% Display the point cloud along with the labels
ax = pcshow(ptCloud.Location,labels);
axis on

labels = unique(labels);
cmap = [[255 255 255];  % White Color (Unlabeled Points)
        [255  0   0];   % Red Color (Dangerous Points)
        [150  75  0];   % Brown Color (Unsuitable Points)
        [255  255 0];   % Yellow Color (Risky Points)
        [0    255 0]];  % Green Color (Suitable Points)
cmap = cmap./255;
cmap = cmap(labels,:);
colormap(ax,cmap)

% Add colorbar to current figure
c = colorbar(ax);
c.Color = "w";
numClasses = numel(labels);

% Center tick labels and use class names for tick marks
c.Ticks = linspace(c.Ticks(1)+(c.Ticks(end)-c.Ticks(1))/(2*numClasses), ...
    c.Ticks(end)-(c.Ticks(end)-c.Ticks(1))/(2*numClasses), ...
    numClasses);
c.TickLabels = classNames(labels);

% Remove tick mark
c.TickLength = 0;
end

The helperNeighborIndices helper function returns the nearest neighbor indices within the specified radius of each point. The function also labels a point as unsuitable when the number of neighbors within the specified radius for the point is less than the specified minimum number of points, or when the point has dangerous points in its neighborhood.

function [curLabels,curNeighborInd] = helperNeighborIndices( ...
    innerBlockPtCloud,outerBlockPtCloud,curLabels, ...
    curInnerBlockLabelsToOuterBlockLabels,minPoints,radius)
curNeighborInd = cell(innerBlockPtCloud.Count,1);
for i = 1:innerBlockPtCloud.Count
    % Compute nearest neighbor for only the unlabeled points
    if isnan(curLabels(curInnerBlockLabelsToOuterBlockLabels(i)))

        % Find the nearest neighbors within radius for each point
        [indices,~] = findNeighborsInRadius(outerBlockPtCloud, ...
            innerBlockPtCloud.Location(i,:),radius);

        % If the number of neighbors is less than the minimum neighbors,
        % label the point as unsuitable
        if numel(indices) < minPoints
            curLabels(curInnerBlockLabelsToOuterBlockLabels(i)) = 2;
            continue
        end

        % If the point is unlabeled and there are dangerous points present in
        % the vicinity of the point, label the point as unsuitable
        if any(curLabels(indices) == 1) && ...
                isnan(curLabels(curInnerBlockLabelsToOuterBlockLabels(i)))
            curLabels(curInnerBlockLabelsToOuterBlockLabels(i)) = 2;
            continue
        end

        curNeighborInd{i} = indices;
    end
end
end

The helperLabelRiskyPts helper function labels the points in the vicinity of the unsuitable points as risky. Additionally, it evaluates a point and its neighbors based on specified attributes. The function labels a point as risky if any of its attribute values is greater than the specified threshold.

function [curLabels,innerBlockIndicesCount, ...
    totalNeighborInd] = helperLabelRiskyPts(innerBlockIndicesCount, ...
    totalNeighborInd,innerBlockPtCloud,outerBlockPtCloud,curLabels, ...
    curInnerBlockLabelsToOuterBlockLabels, ...
    verticalVarianceThreshold,slopeThreshold,residualThreshold, ...
    reliefThreshold)
for i = 1:innerBlockPtCloud.Count
    indices = totalNeighborInd{innerBlockIndicesCount};
    innerBlockIndicesCount = innerBlockIndicesCount + 1;

    if ~isempty(indices) && ...
            isnan(curLabels(curInnerBlockLabelsToOuterBlockLabels(i)))

        % If the point has neighbors, the point is unlabeled. and there
        % are unsuitable points in the vicinity of the point, label the
        % point as risky
        if any(curLabels(indices) == 2)
            curLabels(curInnerBlockLabelsToOuterBlockLabels(i)) = 3;
            totalNeighborInd{innerBlockIndicesCount-1} = [];
            continue
        end

        % Evaluate the point for the vertical variance attribute
        verticalVariance = var(outerBlockPtCloud.Location(indices,3));
        if verticalVariance > verticalVarianceThreshold
            curLabels(curInnerBlockLabelsToOuterBlockLabels(i)) = 3;
            totalNeighborInd{innerBlockIndicesCount-1} = [];
            continue;
        end

        % Evaluate the point for the relief attribute
        relief = max(outerBlockPtCloud.Location(indices,3)) ...
            -min(outerBlockPtCloud.Location(indices,3));
        if relief > reliefThreshold
            curLabels(curInnerBlockLabelsToOuterBlockLabels(i)) = 3;
            totalNeighborInd{innerBlockIndicesCount-1} = [];
            continue
        end

        % Perform the plane fitting operation on the point and its
        % neighbors
        [model,~,outlierIndices] = pcfitplane(pointCloud( ...
            outerBlockPtCloud.Location(indices,:)), ...
            0.3,[0,0,1]);

        % Evaluate the point for the slope attribute
        slope = acosd(model.Normal(3));
        if slope > slopeThreshold
            curLabels(curInnerBlockLabelsToOuterBlockLabels(i)) = 3;
            totalNeighborInd{innerBlockIndicesCount-1} = [];
            continue
        end

        % Evaluate the point for the residual attribute
        residual = rms(abs(outerBlockPtCloud.Location(outlierIndices,:) ...
            *model.Parameters(1:3)' + model.Parameters(4)));
        if residual > residualThreshold
            curLabels(curInnerBlockLabelsToOuterBlockLabels(i)) = 3;
            totalNeighborInd{innerBlockIndicesCount-1} = [];
            continue
        end
    end
end
end

References

[1] OpenTopography. “Tuscaloosa, AL: Seasonal Inundation Dynamics And Invertebrate Communities,” 2011. https://doi.org/10.5069/G9SF2T3K.

[2] Yan, Lu, Juntong Qi, Mingming Wang, Chong Wu, and Ji Xin. "A Safe Landing Site Selection Method of UAVs Based on LiDAR Point Clouds." In 2020 39th Chinese Control Conference (CCC), 6497–6502. Shenyang, China: IEEE, 2020.https://doi.org/10.23919/CCC50068.2020.9189499.