Calculating Normalized Heights¶
Author: | Bradley Chambers |
---|---|
Contact: | brad.chambers@gmail.com |
Date: | 11/11/2015 |
This tutorial will describe the creation of a new filter for calculating
normalized heights, filters.height
.
Note
filters.height
required PCL and has since been replaced by
filters.hag, which is a native PDAL filter. We leave this tutorial as
an example of how to create a filter, and of how to work with the PCL plugin,
but the filter in reference is no longer available.
Introduction¶
Normalized heights are a commonly used attribute of point cloud data. This can also be referred to as height above ground (HAG) or above ground level (AGL) heights. In the end, it is simply a measure of a point’s relative height as opposed to its raw elevation value.
The process of computing normalized heights is straightforward. First, we must have an estimate of the underlying terrain model. With this we can compute the difference between each point’s elevation and the elevation of the terrain model at the same XY coordinate. The quality of the normalized heights will be a function of the quality of the terrain model, which of course depends on the quality of the ground segmentation approach and any interpolation that is required to arrive at the terrain elevation for a given XY coordinate.
We will use a nearest neighbor interpolation scheme to estimate terrain elevations. While this may not be the most accurate approach, it is available in PDAL today, and we hope it will inspire you to implement your own methods!
Approach¶
For the height filter, we only assume that our input point cloud has an already
existing Classification
dimension with some subset of points marked as
ground (Classification=2
). This could, for example, be generated by
filters.pmf (see Identifying ground returns using ProgressiveMorphologicalFilter segmentation), but you can use whichever method
you choose, as long as the ground returns are marked.
Note
We expect ground returns to have the classification value of 2 in keeping with the ASPRS Standard LIDAR Point Classes (see http://www.asprs.org/a/society/committees/standards/LAS_1_4_r13.pdf).
For the full source, please see HeightFilter.cpp
in the
plugins/pcl/filters
folder. Below, we dissect the key aspects of the
algorithm and its implementation.
The bulk of our processing is actually taking place within PCL. For convenience, we’ve defined:
typedef pcl::PointCloud<pcl::PointXYZ> Cloud;
Our first step is to convert the filter’s incoming PDAL PointView
to a PCL
PointCloud
, which requires that we first calculate our bounds so that we
can subtract our XYZ offsets in the conversion step.
BOX3D bounds;
view.calculateBounds(bounds);
Cloud::Ptr cloud_in(new Cloud);
pclsupport::PDALtoPCD(std::make_shared<PointView>(view), *cloud_in, bounds);
Next, we will create two vectors of indices - one for ground returns, one for non-ground returns - and make our first pass through the point cloud to populate these.
pcl::PointIndices::Ptr ground(new pcl::PointIndices());
ground->indices.reserve(view.size());
std::vector<PointId> nonground;
nonground.reserve(view.size());
for (PointId id = 0; id < view.size(); ++id)
{
double c = view.getFieldAs<double>(Dimension::Id::Classification, id);
if (c == 2)
ground->indices.push_back(id);
else
nonground.push_back(id);
}
With our ground indices identified, we can use PCL to extract the ground
returns into a new PointCloud
.
pcl::ExtractIndices<pcl::PointXYZ> extract;
extract.setInputCloud(cloud_in);
extract.setIndices(ground);
Cloud::Ptr cloud_ground(new Cloud);
extract.setNegative(false);
extract.filter(*cloud_ground);
We repeat the extraction now, flipping setNegative
from false to true to
extract the non-ground returns into a new PointCloud
.
Cloud::Ptr cloud_nonground(new Cloud);
extract.setNegative(true);
extract.filter(*cloud_nonground);
To compute the normalized height, we wish to find the nearest ground point for each non-ground point. Here, we achieve this by using a nearest neighbor interpolation scheme. One may prefer to use a more sophisticated interpolation scheme, but that is beyond the scope of this tutorial. We begin by defining model coefficients that will allow us to project the ground and non-ground clouds into the XY plane.
pcl::ModelCoefficients::Ptr coefficients(new pcl::ModelCoefficients());
coefficients->values.resize(4);
coefficients->values[0] = coefficients->values[1] = 0;
coefficients->values[2] = 1.0;
coefficients->values[3] = 0;
We can now project the ground points
pcl::ProjectInliers<pcl::PointXYZ> proj;
proj.setModelType(pcl::SACMODEL_PLANE);
Cloud::Ptr cloud_ground_projected(new Cloud);
proj.setInputCloud(cloud_ground);
proj.setModelCoefficients(coefficients);
proj.filter(*cloud_ground_projected);
followed by the non-ground points
Cloud::Ptr cloud_nonground_projected(new Cloud);
proj.setInputCloud(cloud_nonground);
proj.setModelCoefficients(coefficients);
proj.filter(*cloud_nonground_projected);
Next, we create a KdTree to accelerate our nearest neighbor search. The tree is composed of only ground returns, as our non-ground returns will serve as query points for the nearest neighbor search.
pcl::search::KdTree<pcl::PointXYZ>::Ptr ground_tree;
ground_tree.reset(new pcl::search::KdTree<pcl::PointXYZ> (false));
ground_tree->setInputCloud(cloud_ground_projected);
We iterate over each of our projected non-ground points, searching for our
nearest neighbor in the ground points. Using the indices of each the query
(non-ground) and nearest neighbor (ground), we can retrieve the Z dimension
from the input cloud, compute the height, and set this field in our original
PointView
.
for (int i = 0; i < cloud_nonground_projected->size(); ++i)
{
pcl::PointXYZ nonground_query = cloud_nonground_projected->points[i];
std::vector<int> neighbors(1);
std::vector<float> sqr_distances(1);
ground_tree->nearestKSearch(nonground_query, 1, neighbors, sqr_distances);
double nonground_Z = view.getFieldAs<double>(Dimension::Id::Z, nonground[i]);
double ground_Z = view.getFieldAs<double>(Dimension::Id::Z, ground->indices[neighbors[0]]);
double height = nonground_Z - ground_Z;
view.setField(m_heightDim, nonground[i], height);
}
The only thing left is to set the height field to 0.0 for each of the ground points.
for (auto const& ground_idx : ground->indices)
view.setField(m_heightDim, ground_idx, 0.0);
Example #1¶
Using the autzen dataset (here shown colored by elevation)

we run the following PDAL CLI command
$ pdal translate autzen.laz autzen-height.bpf height \
--writers.bpf.output_dims="X,Y,Z,HeightAboveGround"
The result, when colored by the normalized height instead of elevation is

Example #2¶
If you’d like to overwrite your Z values, follow the height filter with filters.ferry.
$ pdal translate input.laz output-height-as-Z.bpf height ferry \
--writers.bpf.output_dims="X,Y,Z" --filters.ferry.dimensions="HeightAboveGround=Z"
Example #3¶
If you don’t yet have points classified as ground, start with filters.pmf.
$ pdal translate input.laz output-ground-height.bpf ground height \
--writers.bpf.output_dims="X,Y,Z,HeightAboveGround"