This blog references Keqi Zhang's article "A Progressive Morphological Filter for Removing Nonground Measurements From Airborne LIDAR Data"to remove common features such as large buildings and trees.
Inconvenient peeps can be here:Resource Link Download。
1. Introduction
Airborne LiDAR can acquire high-precision terrain measurements over large areas quickly and at low cost. In order to obtain high-precision DTM/DEM, it is necessary to distinguish between ground points (returned directly from the ground) and non-ground points (buildings, vehicles, vegetation) in the measurement points. Many scholars have adopted various methods to perform "point cloud ground point filtering". (These are described in this blog, so I won't go into them again).
2. Morphological Filters
2.1 Expansion/corrosion
Expansion/erosion are two of the basic operations, in general terms expansion/erosion expands/decreases the size of the feature and this combination is used as an open/close operation. For a LiDAR measurement point p(x, y, z), the expansion operation corresponding to the elevation z value at (x, y) can be defined as:
Where (xp, yp, zp) represents the neighboring points of point p, and w is the operation window (can be a one-dimensional "line" or two-dimensional "rectangle/circle/other shapes, etc."). After the expansion operation is completed, the nearest neighbor of p with the largest elevation value in the window w will be output.
Similarly, the erosion operation is to find the nearest neighbor point with the lowest elevation value within a window w of p points. This can be defined by the following equation:
After understanding the expansion/corrosion of these two basic operations, you can form an open/closed operation by simply combining them, in which the open operation is to corrode and then expand the operation, while the closed operation is to expand and then corrode. LiDAR data processing in the use of "open" operation, the processing effect is shown in the figure below:
It can be seen from the figure that the "dotted line" is the surface created by the "erosion" operation, which eliminates the "tree" points, but leaves the large building incomplete. The "solid line" is the surface created by "expanding" the result of the "erosion" operation, and it can be seen that the shape of the large building has been restored. Based on this, we can see that the "open operation" has the ability to remove small features from the ground and retain large features, which is very important for subsequent processing.
2.2 Morphological filtering
The above "open operation" only removes small features and retains large features, but does not remove all non-ground points, and it is not possible to realize the extraction of complex surfaces by only one "open operation". Therefore, how to utilize the ability of the "open operation" for the next step of extraction is the key to further improvement.
Kilian et al. proposed that a "minimum point" can be found in each "window" of the "open operation" result, and then all other points in the window that fall within one elevation of the "minimum point" are considered to be ground points. Other points in the window that fall within an elevation range of the "nadir" are then considered to be ground points. This elevation range is usually defined according to the scanning accuracy of the airborne LiDAR system, and is normally 20-30 cm.
There are two aspects of this method that are very important to the goodness of the final result:
1. Filter window size, if the window size is set small, it can retain the details of the ground undulation very well, but can only remove small features such as cars, trees, etc., and is less effective in removing buildings (roofs are usually considered as the ground). On the other hand, if the window size is set to a larger size, some "ground points" will be over-removed. For example, if a road is adjacent to a small ditch, and the window size is larger than the width of the road, the road may be considered as a non-ground point (because the point in the small ditch has a lower elevation and is considered as the lowest point within the window, while the road is higher and is judged as a non-ground point). are judged to be non-ground points). In addition, some localized hills and dunes will most likely be "cut out".
2. distribution of buildings and trees in specific/localized areas.
Note: Ideally, we could set up a "window" that is small enough to retain ground details. At the same time, it could be large enough to remove buildings, cars, trees, and other features. However, this ideal situation does not exist in the actual dataset.
To solve this problem, Kilian et al. went on to propose that the filtering could be performed multiple times by varying the window size. Each point is given a weight related to the window size, and the larger the window size, the higher the weight of the point. This method gives better results, but does not distinguish ground points from non-ground points in terms of "point level". (The "point level" distinction between ground and non-ground points can be interpolated for better DEM/DTM generation.)
3. Progressive Morphological Filters
From the above analysis in Section 2.2, it is clear that it is difficult for conventional morphological filtering to detect different features with various size variations through a fixed-size window. This problem can be solved by gradually changing the window size.
As shown in the figure below, a window of size l1 is first used to open the original data. The "dotted line" in the figure shows that trees and other features with dimensions smaller than l1 are removed, and terrain features with dimensions smaller than l1 are "excised" (the elevation of the top of the hill is replaced by the smallest value in l1), but buildings with dimensions larger than l1 are retained. Then, the next iteration is performed, the window size is changed to l2, and the result of the previous process is processed by the open operation. The result can be seen from the "solid line", l2 is larger than the size of the building, so the building is also removed, but at the same time, the top of the hill is "excised" to a larger extent. But at the same time, the top of the hill is "removed" to a larger extent.
Something to keep in mind: The problem of removing features of different sizes was solved by gradually increasing the size of the window, but it introduced the problem of "cutting off" parts of terrain features smaller than the window size, such as the tops of "hills".
This emerging problem can be solved by introducing a height difference threshold. The difference in elevation between the roof of a building and a point on the ground is an "abrupt change", whereas the ground elevation is a "gradual change". Therefore, the obvious difference in elevation change between the two can help us make the distinction. Assuming that dhp,1 represents the elevation difference between the original LiDAR measurements and the first iteration of the surface at any given p-point, and dhT,1 represents the elevation difference threshold, then if dhp,1 ≤ dhT,1, point p is considered to be a ground point, and vice versa, if dhp,1 > dhT,1, point p is considered to be a non-ground point. Thereafter, let dhmax(t),1 be the maximum value of the difference between the initial ground point and the filtered surface in the current iteration, then all measurements are retained if dhT,1 > dhmax(t),1 is selected.
Assuming in the second iteration that the maximum elevation difference between the last filtered surface and the current filtered surface is dhmax(t),2, then if dhmax(t),2 < dhT,2, then any ground point with an elevation difference value in the range of dhmax(t),2 will be retained. Similarly, assuming that the building elevation difference between the last iteration and the current iteration is minimized to dhmin(b),2 (which is usually approximated as the height of the building), the building is removed if dhmin(b),2 > dhT,2.
Usually dhT,k is set as the shortest elevation value of the building in the kth iteration of the study area. Using dhT,k as a threshold, for any point p in the kth iteration if dhp,k < dhT,k it is labeled as a ground point and vice versa as a non-ground point. In this way, buildings (trees) of different sizes can be progressively removed as the size of the iteration window increases.
In summary, the detailed flow of Progressive Morphological Filters is shown below:
We can summarize the following four steps for the above diagram:
- Load the irregular point cloud, divide it into regular grids, select the lowest point of elevation in each grid (or interpolate according to the nearest neighbor if there is no point in the grid), and construct the minimum elevation surface.
- The first iteration is performed using the input initial filter window size, and the minimum elevation surface acquired in 1) as parameters for the first iteration. Subsequently, in subsequent iterations, the filtered surface obtained in the previous iteration and the filter window size calculated in 3) are used as inputs. The output of each iteration has two parts: a) a smoother surface obtained after morphological filtering and b) non-ground points detected based on different thresholds.
- Calculate a new filter window size and a different elevation interpolation threshold. Repeat steps 2) and 3) until the window size is larger than the preset maximum window.
- Generation of DEM/DTM based on the removal of non-ground measurements.
Note that the "open operations" in each iteration actually act on the points in the mesh delineated in step 1), so Progressive Morphological
Filters are "point level" filters for LiDAR measurements.
3.1 Parameter calculation (window size/elevation difference threshold)
In step 3) above we have to vary the values of two parameters, window size wk and elevation difference threshold dhT,k, for the next iteration, so how are these two values calculated?
3.1.1 Window size
First is the window size wk is calculated in two ways, the first being:
where k is the number of iterations, b is the initial window size (input by the user), and finally +1 is to ensure that wk is an odd number and the window is symmetric. However, if a study area has very large features, this increasing window size too slowly will take more time. Therefore, a second approach can be used to change the window size by exponential growth, calculated as in the following equation:
Again, where k is the number of iterations and b is the initial window size (input by the user), this approach grows much faster than the first approach.
3.1.2 Elevation difference thresholds
The elevation difference threshold is inextricably linked to the topographic slope of the study area and can therefore be calculated by the following equation:
where dh0 is the initial elevation difference threshold, s is the slope, c is the grid size, and dhmax is the maximum elevation difference threshold.
In urban areas, trees and cars are much smaller in size relative to buildings, so buildings are usually the last to be filtered out, and the maximum elevation difference threshold dhmax can be set to a fixed value (e.g., the height of the shortest building). In mountainous areas, on the other hand, the main non-ground points are vegetation, so there is no need to set a fixed maximum elevation difference threshold dhmax, and so dhmax is usually set to the maximum elevation difference within the survey area.
In addition, the slope s is calculated by the maximum elevation difference dhmax(t),k of the kth iteration, and the window size wk, as shown in the following equation:
3.2 Parameter input/output
3.2.1 Parameter input
- Raw LiDAR point cloud data with each point represented by (x, y, z)
- Division of grid sizec Parameterb (calculation of window size)
- Maximum window size (to determine whether to stop iterating)
- Terrain slope s
- Initial elevation difference threshold dh0
- Maximum elevation difference dhmax
3.2.1 Parameter output
- ground point
- non-terrestrial point
3.3 PCL Official Sample Code
#include <iostream>
#include <pcl/io/pcd_io.h>
#include <pcl/point_types.h>
#include <pcl/filters/extract_indices.h>
#include <pcl/segmentation/progressive_morphological_filter.h>
int
main (int argc, char** argv)
{
pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);
pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_filtered (new pcl::PointCloud<pcl::PointXYZ>);
pcl::PointIndicesPtr ground (new pcl::PointIndices);
// Fill in the cloud data
pcl::PCDReader reader;
// Replace the path below with the path where you saved your file
<pcl::PointXYZ> ("", *cloud);
std::cerr << "Cloud before filtering: " << std::endl;
std::cerr << *cloud << std::endl;
// Create the filtering object
pcl::ProgressiveMorphologicalFilter<pcl::PointXYZ> pmf;
(cloud);
(20);
(1.0f);
(0.5f);
(3.0f);
(ground->indices);
// Create the filtering object
pcl::ExtractIndices<pcl::PointXYZ> extract;
(cloud);
(ground);
(*cloud_filtered);
std::cerr << "Ground cloud after filtering: " << std::endl;
std::cerr << *cloud_filtered << std::endl;
pcl::PCDWriter writer;
<pcl::PointXYZ> ("samp11-utm_ground.pcd", *cloud_filtered, false);
// Extract non-ground returns
(true);
(*cloud_filtered);
std::cerr << "Object cloud after filtering: " << std::endl;
std::cerr << *cloud_filtered << std::endl;
<pcl::PointXYZ> ("samp11-utm_object.pcd", *cloud_filtered, false);
return (0);
}
To this point this article on python point cloud ground point filter (Progressive Morphological Filter) algorithm (PCL library) is introduced to this article, more related python point cloud ground point filter PCL library content, please search for my previous articles or continue to browse the following related articles I hope that you will support me in the future!