# Estimate Stockpile Volume from Aerial Lidar Data

This example shows how to estimate the volume of a stockpile from aerial point cloud data.

*Stockpile* refers to a large supply of materials, such as metals, chemicals, or other inventory, usually held in reserve at a dedicated storage facility. Stockpiles need bulk machinery, such as trucks and bulldozers, for maintenance.

Stockpile measurement techniques help in estimating the weight, volume, and size of a stockpile. Accurate stockpile measurement helps companies in the mining, construction, and shipping industries make profitable business decisions, which include:

Efficient inventory management

Using stockpile data for strategic planning

Protection of personnel, material, and machinery

Systematic downtime and maintenance

The example show how to compute the stockpile volume from an aerial point cloud using these steps.

Load point cloud data.

Extract stockpile points from the point cloud using ground segmentation and plane-fitting.

Transform the stockpile plane to align it with the

*z-*axis.Convert the transformed stockpile point cloud into a surface mesh, and compute the mesh volume.

### Load and Visualize Data

Load the point cloud data into the workspace using the `pcread`

function and visualize the point cloud.

% Read the point cloud ptCloud = pcread("stockpile.pcd"); % Visualize the point cloud figure pcshow(ptCloud) title("Stockpile Point Cloud")

### Segment Ground Plane from Point Cloud

Use these steps to segment the ground plane from the input point cloud.

Define a boundary around the point cloud to identify the ground plane.

Fit a plane to the boundary points.

Identify the boundary points that form the outer edge of the ground surface of the stockpile.

```
% Define a boundary around the point cloud
boundaryIndices = boundary(double(ptCloud.Location(:,1)),double(ptCloud.Location(:,2)));
boundaryPtCloud = pointCloud([ptCloud.Location(boundaryIndices,:)]);
```

Fit a plane to these boundary points by using the `pcfitplane`

function.

% Specify the maximum distance from an inlier point to the plane maxDistance = 1.0; % Fit a plane to the point cloud with boundary points [groundPlane,inlierIndices] = pcfitplane(boundaryPtCloud,maxDistance); groundPlanePtCloud = select(boundaryPtCloud,inlierIndices); % Visualize the ground points figure pcshow(groundPlanePtCloud) title("Ground Plane of Stockpile")

### Transform Point Cloud to XY-Plane

To accurately compute the volume, transform the input point cloud to align it with the *xy-*plane.

First, estimate the rigid transformation between the ground plane of the point cloud and the *xy-*plane by using the `normalRotation`

function.

```
% Estimate the rigid transformation
referenceVector = [0 0 1];
tform = normalRotation(groundPlane,referenceVector);
```

Next, transform the point cloud to align it with the *z-*axis by using the `pctransform`

function.

```
% Transform the point cloud
stockpilePtCloud = pctransform(ptCloud,tform);
```

After the transformation, some ground points of the stockpile point cloud might lie above or below the *xy-*plane. Apply another transformation to translate the outlier points such that all ground points have a z value of `0`

.

**Note: **You must adjust the translation parameters based on the input point cloud.

% Translation the outlier points if abs(stockpilePtCloud.ZLimits(1)) > 0.1 angles = [0 0 0]; translation = [0 0 -stockpilePtCloud.ZLimits(1)]; tform = rigidtform3d(angles,translation); stockpilePtCloud = pctransform(stockpilePtCloud,tform); end % Visualize the transformed stockpile point cloud figure pcshow(stockpilePtCloud) title("Transformed Stockpile Point Cloud")

### Convert Stockpile Point Cloud to Surface Mesh

Estimate connections between the surface points of the stockpile point cloud by using the Delaunay triangulation method.

```
% Estimate the connections between the vertices
DT = delaunayTriangulation(double(stockpilePtCloud.Location(:,1)),double(stockpilePtCloud.Location(:,2)));
```

Use the vertices to generate a surface mesh by using the `surfaceMesh`

function, and visualize the generated mesh.

% Create a surface mesh around the stockpile point cloud mesh = surfaceMesh(double(stockpilePtCloud.Location),DT.ConnectivityList); % Visualize the surface mesh surfaceMeshShow(mesh)

### Estimate Volume of Stockpile

Estimate the volume of the stockpile mesh.

% Estimate the volume of the stockpile mesh volume = 0; for i = 1:size(mesh.Faces,1) tri = mesh.Faces(i,:); x = mesh.Vertices(tri(1),:); y = mesh.Vertices(tri(2),:); z = mesh.Vertices(tri(3),:); partialVol = (x(3)+y(3)+z(3))*(x(1)*y(2)-y(1)*x(2)+y(1)*z(2)-z(1)*y(2)+z(1)*x(2)-x(1)*z(2))/6; volume = volume + partialVol; end disp("Estimated Volume of Stockpile = " + volume + " cubic metres")

Estimated Volume of Stockpile = 10227.4239 cubic metres

### Tips

For accurate volume measurement:

As a preprocessing step to this workflow, remove the outliers and nonground objects from the input point cloud by using functions such as

`segmentGroundSMRF`

.A high-density point cloud can improve the accuracy of volume estimation, especially for irregular or complex stockpiles.