(十五)pcl-filter篇-Progressive Morp

参考文章：
博客：https://blog.csdn.net/AmbitiousRuralDog/article/details/80268920
论文：http://users.cis.fiu.edu/~chens/PDF/TGRS.pdf
官方开运算和闭运算：https://www.sciencedirect.com/topics/computer-science/morphological-operator
中文开运算和闭运算：https://www.cnblogs.com/pingwen/p/12240890.html
源代码：https://github.com/PointCloudLibrary/pcl/blob/master/segmentation/include/pcl/segmentation/impl/progressive_morphological_filter.hpp
点云语义分割：https://arxiv.org/pdf/1911.11236.pdf

论文解读：《A Progressive Morphological Filter for Removing
Nonground Measurements From
Airborne LIDAR Data》

Abstract—Recent advances in airborne light detection and ranging (LIDAR) technology allow rapid and inexpensive measurements of topography over large areas. This technology is becoming a primary method for generating high-resolution digital terrain models (DTMs) that are essential to numerous applications such as flood modeling and landslide prediction. Airborne LIDAR systems usually return a three-dimensional cloud of point measurements from reflective objects scanned by the laser beneath the flight path. In order to generate a DTM, measurements from nonground features such as buildings, vehicles, and vegetation have to be classified and removed. In this paper, a progressive morphological filter was developed to detect nonground LIDAR measurements. By gradually increasing the window size of the filter and using elevation difference thresholds, the measurements of vehicles, vegetation, and buildings are removed, while ground data are preserved. Datasets from mountainous and flat urbanized areas were selected to test the progressive morphological filter. The results show that the filter can remove most of the nonground points effectively.
摘要:机载光探测和测距(LIDAR)技术的最新进展允许对大面积地形进行快速而廉价的测量。 该技术正在成为生成高分辨率数字地形模型(dms)的主要方法，而dms对于洪水建模和滑坡预测等众多应用都是必不可少的。 机载激光雷达系统通常返回由飞行路径下的激光扫描的反射物体的三维点测量云。 为了生成DTM，必须对建筑物、车辆和植被等非地面特征进行分类和移除。 本文提出了一种渐进形态滤波器来检测非地面激光雷达测量数据。 通过逐渐增大滤波器的窗口大小和使用高程差阈值，去除车辆、植被和建筑物的测量值，保留地面数据。 选取山区和平坦的城市地区的数据集，对渐进形态滤波器进行测试。 结果表明，该滤波器能有效地去除大部分非接地点。

Index Terms—Airborne laser altimetry, digital terrain model (DTM), light detection and ranging (LIDAR) data filtering.
关键字：机载激光测高，数字地形模型(DTM)，光探测和测距(LIDAR)数据滤波。

INTRODUCTION
HIGH-RESOLUTION digital terrain models (DTMs) are essential for many geographic information system(GIS)-related analysis and visualization. The airborne light detection and ranging (LIDAR) technology is revolutionizing our way to acquire a high-resolution DTM by allowing rapid and inexpensive measurements of topography over a large area. Airborne LIDAR systems usually obtain measurements for the horizontal coordinates (x,y) and elevation (z) of the reflective objects scanned by the laser beneath the flight path. These measurements generate a three-dimensional cloud of points with irregular spacing. The laser-scanned objects include buildings, vehicles, vegetation (canopy and understory), and “bare ground.” To generate a DTM, measurements from ground and nonground features have to be identified and classified.Removing nonground points from LIDAR datasets has proven to be a challenging task.
介绍
高分辨率数字地形模型(DTMs)是许多地理信息系统(GIS)相关分析和可视化的基础。 机载光探测和测距(LIDAR)技术通过允许对大面积地形进行快速、廉价的测量，正在彻底改变我们获取高分辨率DTM的方式。 机载激光雷达系统通常获取飞行路径下激光扫描的反射物体的水平坐标(x,y)和高度(z)的测量值。 这些测量产生了不规则间距的三维点云。 激光扫描的对象包括建筑物、车辆、植被(冠层和林下植被)和“裸地”。 为了生成DTM，必须对来自地面和非地面特征的测量数据进行识别和分类。 从激光雷达数据集中删除非地面点被证明是一项具有挑战性的任务。
Kraus and Pfeifer [1], [2] utilized linear least squares interpolation iteratively to remove tree measurements and generate DTMs in forest areas. This method was extended later to filter buildings and trees in urban areas by Pfeifer et al. [3]. The iterative linear interpolation method removes a low-degree polynomial trend surface from the original elevation data to produce a set of reduced elevation values. This method requires that the reduced elevation follows a random process of ergodic property. However, this property is hard to satisfy in urban areas where significant anthropogenic modification of natural terrain occurs. Therefore, the iterative linear interpolation is not guaranteed to converge when being applied to LIDAR measurements for these areas.
Kraus和Pfeifer[1]、[2]采用线性最小二乘插值迭代去除树木测量值，生成森林区域的dtm。 该方法后来被Pfeifer et al.[3]推广到城市地区的建筑物和树木的过滤。 迭代线性插值方法从原始高程数据中去除一个低次多项式趋势面，生成一组降低的高程值。 这种方法要求降低的高程遵循一个随机的遍历性过程。 然而，在自然地形发生重大人为改变的城市地区，这种特性很难满足。 因此，在对这些区域进行激光雷达测量时，不能保证迭代线性插值的收敛性。
Vosselman [4] proposed a slope-based filter that identifies ground data by comparing slopes between a LIDAR point and its neighbors. A point is classified as a ground measurement if the maximum value of slopes between this point and any other point within a given circle is less than a predefined threshold. The lower the threshold slope, the more objects will be removed. The threshold slope for a certain area is either constant or a function of distance. A reasonable threshold slope can be obtained by using prior knowledge about terrain in the study area.
Vosselman[4]提出了一种基于坡度的滤波器，通过比较LIDAR点与其相邻点之间的坡度来识别地面数据。 如果一个点与给定圆内任何其他点之间的斜率最大值小于预定义的阈值，则该点被归类为地面测量。 阈值斜率越低，被移除的对象越多。 某一区域的阈值斜率要么是常数，要么是距离的函数。 利用研究区地形的先验知识，可以得到合理的阈值坡度。
There are two basic errors in classifying LIDAR measurements by virtually any filtering method. One is commission error that classifies nonground points as ground measurements[5]. The other is omission error that removes ground points mistakenly. The critical step in slope-based filtering is to determine an optimum threshold so that omission and commission errors can be minimized. Determining a slope threshold in terms of terrain information in the analyzed area is somewhat subjective. Vosselman [4] demonstrated that good results could be obtained by using threshold slopes from training datasets. However, the training datasets have to include all types of ground measurements in a study area to achieve good results, which is not always practical. Other attempts to improve the slope-based filter can be found in [6] and [7].
有两个基本的错误分类激光雷达测量几乎任何滤波方法。 一种是将非接地点分类为地面测量[5]的误差。 另一种是误移接地点的遗漏误差。 斜率滤波的关键步骤是确定最优阈值，使遗漏和调试误差最小。 根据所分析区域的地形信息确定坡度阈值有点主观。 Vosselman[4]证明了利用训练数据集的阈值斜率可以获得较好的结果。 然而，训练数据集必须包括研究区域内所有类型的地面测量数据才能取得良好的效果，这并不总是可行的。 改进基于坡度的滤波器的其他尝试可以在[6]和[7]中找到。
The implicit premise of applying the slope-based filter is that there is a distinct difference between the slope of terrain and that of nonground objects such as trees and buildings. Satisfactory results have been achieved in our experiments by applying the slope-based filter to flat urban areas such as Miami, FL. However, both omission and commission errors were large when this method was applied to vegetated mountain areas with a large slope variation.
应用基于坡度的滤波的隐含前提是，地形的坡度与树木、建筑物等非地物的坡度有明显的区别。 在我们的实验中，将基于坡度的滤波方法应用于较为平坦的城市地区，如迈阿密，FL，取得了令人满意的结果。然而，当该方法应用于植被覆盖的坡度变化较大的山区时，漏检误差和执行误差都较大。
Haugerud and Harding [8] developed an algorithm to filter tree points in forest areas by comparing local curvatures of point measurements. Ground measurements were selected by removing tree vertices iteratively from a triangulated irregular network (TIN) constructed from LIDAR measurements. Alternatively, ground points can be classified by iteratively selecting ground measurements from an original dataset. Axelsson [9] suggested adaptive TIN models to find ground points in urban areas. First, seed ground points within a user-defined grid of a size greater than the largest nonground features are selected to compose an initial ground dataset. Then, one point above each TIN facet is added to the ground dataset every iteration if its parameters are below threshold values. The iteration continues until no points can be added to the ground dataset. The problem with the adaptive TIN method is that different thresholds have to be given for various land cover types.
Haugerud和Harding[8]开发了一种算法，通过比较点测量值的局部曲率来过滤森林区域中的树点。 地面测量是通过从由激光雷达测量构建的不规则三角网络(TIN)中迭代去除树顶点来选择的。 另外，可以通过从原始数据集迭代地选择地面测量值来对地面点进行分类。 Axelsson[9]提出了自适应TIN模型来寻找城市地区的接地点。 首先，选择用户定义网格中大于最大非地面特征的种子地点组成初始地面数据集。 然后，如果TIN facet的参数低于阈值，则在每次迭代时将每个TIN facet上面的一个点添加到地面数据集。 迭代将继续进行，直到不能向地面数据集添加点为止。 自适应TIN方法的问题是，对于不同的土地覆盖类型，必须给出不同的阈值。
Another commonly used algorithm to remove nonground objects is a mathematical morphology filter which is applied to a grayscale image [4], [10], [11]. The elevation of trees, cars,and buildings is usually higher than those of surrounding ground points. If LIDAR points are converted into a regular, grayscale grid image in terms of elevation, then the shapes of buildings, cars, and trees can be identified by the change of gray tone. It is well known that compositions of algebraic set operations based on mathematical morphology can be used to identify objects in a grayscale image [12]. Therefore, mathematical morphology can be used to filter LIDAR data. The main objective of this paper is to develop a progressive morphology filter to enable automatic extraction of ground points from LIDAR measurements with minimal human interaction.
另一种常用的去除非地物的算法是数学形态学滤波器，它应用于灰度图像[4]，[10]，[11]。 树木、汽车和建筑物的标高通常高于周围的地面点。 如果按照高程将LIDAR点转换成规则的灰度网格图像，则可以通过灰度色调的变化来识别建筑物、汽车、树木的形状。 众所周知，基于数学形态学的代数集运算组合可以用于识别灰度图像[12]中的目标。 因此，数学形态学可以用来过滤激光雷达数据。 本文的主要目标是开发一个渐进的形态滤波器，使自动提取地面点从激光雷达测量与最小的人的互动。 （人工干预）
Mathematical morphology composes operations based on set theory to extract features from an image. Two fundamental operations, dilation and erosion, are commonly employed to enlarge (dilate) or reduce (erode) the size of features in binary images. Dilation and erosion operations may be combined to produce opening and closing operations. The concept of erosion and dilation has been extended to multilevel (grayscale) images and corresponds to finding the minimum or maximum of the combinations of pixel values and the kernel function, respectively, within a specified neighborhood of each raster [13]. These concepts can also be extended to the analysis of a continuous surface such as a digital surface model as measured by LIDAR data. For a LIDAR measurement p(x,y,z), the dilationof elevation z at x and y is defined as

(1).png
数学形态学是一种基于集合理论的运算，用于从图像中提取特征。 在二值图像中，通常采用膨胀和侵蚀两种基本运算来放大或减小特征的尺寸。 膨胀和侵蚀作业可以结合起来进行开合作业。 侵蚀和膨胀的概念已扩展到多级(灰度)图像，对应于在每个光栅[13]的指定邻域内分别找到像素值和核函数组合的最小或最大值。 这些概念也可以扩展到连续表面的分析，如激光雷达数据测量的数字表面模型。 对于一个LIDAR测量p(x,y,z)，高度z在x和y处的膨胀定义为 (1).png

where points (xp,yp,zp) represent p's neighbors (coordinates) within a window,w . The window can be a one-dimensional(1-D) line or two-dimensional (2-D) rectangle or other shapes. The dilation output is the maximum elevation value in the neighborhood of p. Erosion is a counterpart of dilation and is defined as

(2).png
其中点(xp,yp,zp)表示窗口中p的邻居(坐标)，w。 窗口可以是一维(1-D)线或二维(2-D)矩形或其他形状。 膨胀输出为p附近的最大高程值。侵蚀是膨胀的对应物，定义为 (2).png
Fig2.png
Fig.2 Process of the progressive morphological filter to identify terrain and building measurements. The dots represent synthetic points based on LIDAR surveys. The first filtered elevation surface (dashed line) is obtained by applying an opening operation with a window size of 15 m (l1 ) to the raw point data. The second filtered elevation surface (solid line) is derived by applying an opening operation with a window size of 21 m (l2) to the first filtered surface.
Fig2.png

图2过程的渐进形态滤波器，以识别地形和建筑测量。 这些点表示基于激光雷达调查的合成点。 第一个过滤的高程面(虚线)是通过对原始点数据应用一个窗口大小为15 m (l1)的开放操作获得的。 第二过滤的高程面(实线)是通过对第一过滤的表面应用一个窗口大小为21米(l2)的开孔操作得到的。

Fig. 2 illustrates the process of a progressive morphological filter. An initial filtered surface is derived by applying an opening operation with a window of length l1 to the raw data. The large nonground features such as buildings are preserved because their sizes are larger than l1, while individual trees of size smaller than l1 are removed. For the terrain, features smaller than l1 are cut off and replaced by the minimum elevation within l1.
图2显示了渐进式形态滤波器的过程。 通过对原始数据应用长度为l1的窗口操作，得到了一个初始滤波曲面。 因为建筑物的大小大于l1，所以保存了大型的非地面特征，而小于l1的树则被去除。 对于地形，将小于l1的特征截断，用l1内的最小标高代替。
In the next iteration, the window size is increased to l2, and another opening operation is applied to the filtered surface, resulting in a further smoothed surface. The building measurements are removed and replaced by the minimum elevation of previous filtered surface within l2, since the size of the building is smaller than the current window size.
在下一次迭代中，窗口大小增加到l2，并对过滤后的表面进行另一个开孔操作，使表面进一步平滑。 由于建筑的尺寸比当前的窗户尺寸小，所以建筑的测量值被l2内先前过滤过的表面的最小标高所取代。

By performing an opening operation to laser-scanned data with a line window that increases in size gradually, the progressive morphological filter can remove buildings and trees at various sizes from a LIDAR dataset. However, the filtering process tends to produce a surface that lies below the terrain measurements, leading to incorrect removal of the measurements at the top of high-relief terrain (points above l2 in Fig. 2). Even in the flat ground areas, the filtered surface is usually lower than the original measurements. Therefore, most point measurements for terrain are removed, and only a filtered surface is available if the opening operation is performed to the LIDAR data directly. This problem can be overcome by introducing an elevation difference threshold based on elevation variations of the terrain, buildings,and trees.
通过对激光扫描数据执行一个开操作，该行窗口的大小逐渐增加，渐进式形态滤波器可以从LIDAR数据集中去除各种大小的建筑和树木。 然而，滤波过程往往会产生一个低于地形测量值的表面，导致高起伏地形顶部(图2中l2以上点)的测量值被不正确地移除。即使在平坦的地面区域，滤波后的表面通常也低于原始测量值。 因此，对于地形的大多数点测量都被删除了，如果直接对LIDAR数据进行打开操作，则只有一个过滤过的表面可用。 该问题可以通过引入基于地形、建筑物和树木的高程变化的高程差阈值来克服。
Each building has a certain size and height. There is an abrupt change in elevation between the roof and base of a building, while the elevation changes of terrain are gradual (Fig. 2). The difference in the elevation variations of buildings and terrain can help the filter to separate the building and ground measurements. Suppose that dhp,1 represents the height difference between an original LIDAR measurement and the filtered surface in the initial iteration at any given point p(Fig. 2), and dhT,1 represents the elevation difference threshold. Point p is classified as a ground measurement if dhp,1 <= dhT,1 and as a nonground measurement if dhp,1 > dhT,1. Let dhmax(t),1 stand for the maximum height difference between the original terrain measurements and the filtered surface (Fig. 2). If a dhT,1 is selected such that the dhmax(t),1 is less than dhT,1, then the LIDAR measurements for terrain will be preserved. In general, dhT,1 will be a function of window size. How to derive this threshold value will be discussed later.
每座建筑都有一定的大小和高度。 建筑物的屋顶与地基之间存在高程突变，而地形的高程变化是渐变的(图2)。建筑物与地形高程变化的差异可以帮助过滤器将建筑物与地面测量分开。 假设dhp,1表示初始迭代中任意给定点p处的原始LIDAR测量值与被过滤表面的高度差(图2)。 2)， dhT,1为高程差阈值。 当dhp,1 <= dhT,1时，p为地面测量值;当dhp,1 > dhT,1时，p为非地面测量值。 设dhmax(t)，1表示原始地形测量值与滤波后表面的最大高差(图2)，如果选择dhT,1，使dhmax(t)，1小于dhT,1，则保留地形的LIDAR测量值。 一般来说，dhT,1是窗口大小的函数。 稍后将讨论如何推导这个阈值。
In the second iteration, suppose that the maximum height difference between the previous and this filtered terrain surface is dhmax(t),2. The ground measurements within dhmax(t),2 will be preserved as long as dhmax(t),2 is smaller than the elevation difference threshold dhT,2 for the current operation. Suppose that the minimum elevation difference for the building between the previous and current filtering operation is represented by dhmin(b),2, which is approximately the height of the building.The building measurements will be removed on condition that dhmin(b),2 is larger than dhT,2.
在第二次迭代中，假设前一个过滤后的地形表面与该地形表面的最大高差为dhmax(t)，2。dhmax(t)，2以内的地面测量值，只要dhmax(t)，2小于当前操作的高程差阈值dhT,2，就会保留。假设上述滤波操作与当前滤波操作的建筑物的最小高差用dhmin(b)，2表示，其近似为建筑物的高度。当dhmin(b)，2大于dhT,2时，建筑测量值将被移除。
Generally, the elevation difference threshold dhT,k is set to be the minimum height value of the building objects in an analyzed area at iteration k. Taking dhT,k as the threshold, for any given point p at k.th opening operation, we mark p as a ground measurement if dhp,k <= dhT,k , and as a nonground measurement otherwise. In this way, the measurements for buildings with various sizes can be identified by gradually increasing the widow sizes and applying an opening operation repeatedly until a window size is greater than the size of the largest building. Since there is also an abrupt elevation change from trees to adjacent ground, the above building filtering procedure can be applied to the removal of tree measurements as well. Note that the filtered surfaces from the opening operation are not utilized to generate the DTM, but used to help classify point measurements together with elevation difference thresholds.
通常,高差阈值dhT, k是将建筑物的最小高度值对象在一个分析区域迭代k。dhT, k为阈值,对于任何给定的点p在k.th开放操作,我们p标记为地面测量如果dhp,k <= dhT,k ,否则nonground测量。这样，可以通过逐渐增大窗窗大小，并反复进行开窗操作，直到窗口大小大于最大建筑的大小，来确定不同大小建筑的测量值。由于树木与相邻地面之间也存在高程突变，因此上述建筑物过滤程序也可用于去除树木测量值。注意，从开口操作中过滤出来的表面并没有被用来生成DTM，而是用来帮助分类点测量值和高程差阈值。
The detailed steps to use the progressive morphological filter to construct the DTMs are shown in Fig. 3 and given as follows.
使用渐进式形态滤波器构造dms的具体步骤如图3所示，如下所示。

Fig.3.png
Step 1: The irregularly spaced (x; y; z) LIDAR measurements are loaded. A regularly spaced minimum surface grid is constructed by selecting the minimum elevation in each grid cell. Point coordinates (x; y; z) are stored in each grid cell. If a cell contains no measurements, it is assigned the value of nearest point measurement.
Step 2: The progressive morphological filter whose major component is an opening operation is applied to the grid surface. At the first iteration, the minimum elevation surface together with an initial filtering window size provide the inputs for the filter. In the following iterations, the filtered surface obtained from the previous iteration and an increased window size from Step 3 are used as input for the filter. The output of this step include a) the further smoothed surface from the morphological filter and b) the detected nonground points based on the elevation difference threshold.
Step 3: The size of the filter window is increased and the elevation difference threshold is calculated. Steps 2 to 3 are repeated until the size of the filter window is greater than a predefined maximum value. This value is usually set to be slightly larger than the maximum building size.
Step 4: The last step is to generate the DTMs based on the dataset after nonground measurements have been removed.
步骤1:不规则间隔(x;y;z)加载激光雷达测量数据。通过选择每个网格单元的最小高程，构造出规则间隔的最小曲面网格。点坐标(x;y;Z)存储在每个网格单元中。如果单元格不包含测量值，则为其分配最近点测量值。
步骤2:将主要成分为开孔操作的渐进式形态滤波器应用于网格表面。在第一次迭代中，最小高程面和初始滤波窗口大小为滤波器提供输入。在接下来的迭代中，从前面的迭代中获得的经过过滤的表面和从步骤3中增加的窗口大小被用作过滤器的输入。这一步的输出包括a)形态滤波器进一步平滑的表面和b)基于高程差阈值检测到的非地面点。
步骤3:增大滤波窗口大小，计算高程差阈值。重复步骤2到3，直到筛选窗口的大小大于预定义的最大值。这个值通常设置为略大于最大建筑尺寸。
步骤4:最后一步是去除非地面测量数据后，根据数据集生成DTMs。
Note that each cell of the minimum surface grid generated in Step 1 contains an original or interpolated LIDAR point with elevation representing the cell value. Any filtering operation performed to the grid is actually applied to points in cells. Therefore, the progressive morphological filter classifies the LIDAR measurements at the point level.
请注意，在步骤1中生成的最小表面网格的每个单元包含一个原始的或插值的激光雷达点，高程代表单元值。对网格执行的任何过滤操作实际上应用于单元格中的点。因此，渐进式形态滤波器在点水平上对激光雷达测量数据进行分类。
IV. PARAMETERS FOR PROGRESSIVE MORPHOLOGICAL FILTER
The selection of the window size and elevation difference threshold is critical to achieve good results when applying the morphological filter. For window size selection, one straightforward choice is to increase the window size linearly by the following formula:
wk = 2kb + 1 (4)
在形态滤波中，窗口大小和高程差阈值的选择是获得良好效果的关键。对于窗口大小的选择，一个直接的选择是通过下面的公式线性增加窗口大小:
wk = 2kb + 1 (4)
where k=1,2,.. , M and b is the initial window. The maximum window size (number of cells) is equal to 2Mb+1 . Taking 2kb+1 as the window size guarantees that the filter window is symmetric around the central point so that the programming of the opening operation is simplified. The advantage of increasing the window size linearly is that gradually changing topographic features are well preserved. However, considerable computing time is needed for an area with large nonground objects. Alternatively, the window size can be increased exponentially to reduce the number of iterations used in the filter.

， M和b为初始窗口。最大窗口大小(单元格数)等于2Mb+1。取2kb+1作为窗口大小，保证了滤波器窗口围绕中心点对称，从而简化了打开操作的编程。线性增加窗口大小的优点是，逐渐变化的地形特征被很好地保存下来。然而，对于一个有大型非地面物体的区域来说，需要相当长的计算时间。或者，窗口大小可以按指数级增加，以减少过滤器中使用的迭代次数。
wk=2b(k次方)+1 （5）
where b is the base of an exponential function; k = 0,1,2,..., M and again 2b(M次方) + 1 is equal to the maximum window size.The elevation difference threshold can be determined based on the slope of topography in the study area. There is a relationship between the maximum elevation difference dhmax(t),k for the terrain, window size wk, and the terrain slope s (Fig. 2) assuming that the slope is constant.

k = 0, 1, 2,…再次,M和2 b (M次方)+ 1等于最大的窗口大小。高差阈值可根据研究区地形坡度确定。假设坡度不变，最大高差dhmax(t)、地形k、窗口大小wk和地形坡度s之间存在关系(图2)。