- Email: [email protected]
An improved detection method based on morphology and profile analysis for bridge extraction from Lidar
Optics and Laser Technology 121 (2020) 105790
Contents lists available at ScienceDirect
Optics and Laser Technology journal homepage: www.elsevier.com/locate/optlastec
Full length article
An improved detection method based on morphology and profile analysis for bridge extraction from Lidar
T
⁎
Lin Lia, Wei Rongb,a, , Fei Sub, Xiaoyu Xingb a b
Collaborative Innovation Centre of Geospatial Technology, Wuhan University, 129 Luoyu Road, Wuhan 430079, China School of Resource and Environmental Sciences, Wuhan University, 129 Luoyu Road, Wuhan 430079, China
H I GH L IG H T S
dimensional discrete points morphological operator, removing vegetation. • Three sets incorporated with the profile analysis method, removing large background objects. • Union-find optimized profile analysis method with the topological characteristics, extracting the bridges. • The • OpenMP to improve the computation efficiency.
A R T I C LE I N FO
A B S T R A C T
Keywords: Bridge extraction Three dimensional morphology Profile analysis OpenMP
Extraction of bridges from light detection and ranging (Lidar) images is a difficult problem with low detection accuracy and detection efficiency because of strong dependence on bridge shapes, the influence of vegetation and the large amount of data. This paper proposes an improved method based on morphology and profile analysis to extract bridges with removing background objects such as vegetation. To remove vegetation effectively, a new morphological three dimensional (3D) discrete points operator is applied which the gridding method is utilized to obtain the neighborhoods of the discrete points. Subsequently, union-find sets is incorporated with the profile analysis method to optimize the process of generating the minimum spanning tree (MST) and determining connected domains, which are used to remove large objects and leave the bridges behind. By combining the optimized profile analysis method with the topological characteristics of bridges, bridges are extracted without dependence on their shapes. Finally, in order to improve the computation efficiency of the Lidar data, the OpenMP is employed. Experimental results show that the proposed method can extract the bridges from Lidar data effectively.
1. Introduction Light detection and ranging (Lidar) is a new airborne remote sensing technology that can be used to quickly and accurately obtain the position and elevation, i.e., the point cloud data, of ground objects. Currently, the main areas in which Lidar is applied include aspects of building extraction and reconstruction [1,2], digital elevation model (DEM) generation [3,4], and residential area recognition based on Lidar and remote sensing images [5]. The bridge not only play an important role in city planning but are also is a crucial military targets in war. Consequently, accurate and rapid detection of bridge has important research value. Most existing bridge detection methods are based on optical remote sensing aerospace images instead of Lidar data. These methods can be classified into
⁎
object-oriented analysis methods [6–8], and knowledge-based methods [9,10]. These methods were developed for the case that the bridges are over water, however, the bridges cannot efficiently extract bridges on land in general. This paper proposes an improved method based on Sithole and Vosselman’s algorithm that overcomes the problems outlined above. The proposed method gradually removes noise points, vegetation, and large objects to extract the bridges from Lidar images accurately. The gridding method is applied for removing the vegetation quickly by a newly designed 3D morphological operator. The profile analysis method is optimized using union-find sets and topology characteristics, through the optimized method, the large background objects can be removed and the bridges are extracted accurately. In order to improve the efficiency of the method further, the parallel processing with
Corresponding author. E-mail addresses: [email protected] (W. Rong), [email protected] (F. Su), [email protected] (X. Xing).
https://doi.org/10.1016/j.optlastec.2019.105790 Received 5 November 2018; Received in revised form 8 May 2019; Accepted 25 August 2019 0030-3992/ © 2019 Elsevier Ltd. All rights reserved.
Optics and Laser Technology 121 (2020) 105790
L. Li, et al.
Fig. 1. Process flow of the proposed method of bridge extraction in this paper.
each other. As a common image processing theory, mathematical morphology is usually used for ground point extraction and point cloud filtering in Lidar data processing [15,16]. In this part of the process, the proposed 3D morphological opening operator can be utilized to remove most of the vegetation and some small objects without affecting large objects. Assuming that A represents a point set from 2D images, for each point a, where a ∈ A, the structure element B is defined as the set of neighborhood points of a. The morphological opening operation is carried out by first using B to corrode A and then expanding the corrosion result with B, as shown in Eq. (1) [17]:
OpenMP is used to obtain the stripe information from various directions in profile analysis. 2. Methods The proposed bridge extraction method consists of four main steps: preprocessing, vegetation removal, large objects removal, and bridge detection as shown in Fig. 1. Following preprocessing with denoising and filtering, a new 3D discrete points morphological operator is applied to remove vegetation. Next, large background objects are removed based on the optimized profile analysis method. Finally, combined with the topology characteristics of bridges, the optimized profile analysis method is used to extract the bridges.
(1)
A°B = (AΘB ) ⊕ B
To apply opening operations on Lidar data, we redefine A as the set of point cloud data of Lidar. For each point a, where a ∈ A, set B is the neighborhood points of a in 3D data space, and a.z represents the height of point a. The corrosion operation on A assigns the height of the lowest point in B to a.z, while the expansion operation on A assigns the height of the highest point in B to a.z. The opening operation on A can be described as:
2.1. Preprocessing Lidar data are obtained from the Lidar system, which comprises a global positioning system (GPS), scanner, inertial measurement unit (IMU), and imaging devices. In general, it uses a laser light to scan the earth and other objects, by calculating the echo time, the equipment can obtain detailed points information about the surface of the ground or the objects, such as height and coordinates. However, the laser light may strike special objects, such as dust in the air, glasses, and metallic shimmering surfaces. In such situations, reflection and refraction of the light will happen. The points of these special data are defined as noise points. Further, owing to the strong penetrability of laser light, the redundant information caused by multiple reflections may exist in Lidar data also. So, it is necessary to filter the noise points and redundant information from the original Lidar data before subsequent operations.
AΘB = \{ a ∈ A|a. z = min (B. z )\} A ⊕ B = \{ a ∈ A|a. z = max (B. z )\} A°B = (AΘB ) ⊕ B
(2)
In general, the neighborhood of one point can be easily obtained in 2D images according to the coordinates. In contrast, in Lidar data, owing to the huge amount of 3D discrete points, the tremendous calculations on the distances between the current point and all the other points significantly reduces the speed of operations. Thus, the gridding method is utilized to rapidly obtain the neighborhood of 3D discrete points. Instead of directly calculating the neighborhood of the current point, by first gridding the Lidar data and ascertaining the neighborhood of the grid containing the current point, the neighborhood of the current point can be determined from the current grid and its neighborhood. The procedure is as follows:
2.1.1. Noise filtering The noise points derived from dangling objects often cause some abnormal data, in which the height differences between the data collector and the ground or the objects are greater than the real value. A simple and efficient way to remove these noise points is to set an appropriate height threshold. For example, if the real height difference is 500 m, we can assign 500 as the threshold, which means data greater than 500 will be discarded.
1. Grid the discrete points (Fig. 2(b)) with grid width twice the average point space, to ensure that all the neighborhood points (usually eight points) are included in the current grid and its neighborhood; here the current grid is defined as the grid containing the current point (point 0). 2. Calculate the distances between all the other points in the nine grids and the current point, and select the eight nearest points using the Top-K sorting algorithm.
2.1.2. Single and last echo data filtering As the penetrability of laser light is very high, when the light reaches an object, especially vegetation, multiple reflections occur. In these reflections, the last echo data usually come from ground. So, selecting the last and single echo data can filtering a certain amount of vegetation and noise points, and also reduce the volume of Lidar data. However, buildings and other large objects do not produce multiple reflections because of their special materials, so this step has no effect on the points representing them. Based on the Lidar standard data format (LAS) [14], single and last echo data filtering can be used to extract those points whose echo numbers are equal to their echo times. After preprocessing, the Lidar data only consist of the points representing the remaining vegetation, large objects, and ground.
This method significantly reduces the search scope and the amount of calculation. Note that the average point space mentioned in step (1) is calculated based on the data format of Lidar points using Eq. (3):
AveragePointSpacing =
2.2. Vegetation removal
(max X − min X )(max Y − min Y ) PointsNum
(3)
where maxX, minX, maxY, and minY are parameters of the data scope, and PointsNum is the overall number of points. Fig. 3 gives an example of vegetation removal with the above 3D morphological opening operations. In the figure, the Lidar data include vegetation, buildings, and roads. The number of points is 98696, and the average point space is 0.75 m. The result in Fig. 3(b) shows that the proposed 3D opening operator removed a reasonable amount of vegetation and some other small objects.
In Lidar data, the points on the same large object have virtually the same height, except for the boundary of the object—especially the area connected with ground, where the points usually have large height differences from the ground. In other words, the points on one large object exist as “patches.” As the ground is usually next to vegetation or some small objects, the points often have greater distances between 2
Optics and Laser Technology 121 (2020) 105790
L. Li, et al.
Fig. 2. Procedure for obtaining the neighborhood of the discrete points using the gridding method: (a) Discrete Lidar Points; (b) Points Gridded; (c) Eight Neighborhood Points (1–8) of Point 0.
constructing Delaunay triangulations and generating the minimum spanning tree (MST). The height threshold and neighborhood threshold, generally set as one-third of average point space and twice the average point space, respectively, are used to delete the edges with high weights in the MST. (3) Overlap the four profiles, which results in some stripes in different directions being connected into a larger connected domain. (4) Determine all the connected domains, each of which represents an object. To accurately indicate the feature of lines and points in the connected domain, six different shapes of lines and points are defined [12], as shown in Fig. 5. For example, if connected region ① only has one adjacent connected domain on the left or right side, and ① is the lower region, then for the shapes of connected region ①, all the lines and points in ① are defined as lowered. Based on the shape of the lines and points in each connected domain [11], set the threshold to determine whether the connected domain belongs to a large object. If it does, then all the points in that connected domain will be directly removed. (5) Iterate steps (1)–(4) between one and four times, depending on the actual components of the point cloud data.
2.3. Large objects removal After filtering the noise points and small objects, the remaining Lidar data mainly consist of the points representing the large objects and the ground. Although existing methods based on the geometric characteristics of Lidar points [3,13,18] perform well in extracting ground points to generate a DEM, they often remove bridges as well when removing large objects. In contrast, the proposed method utilizes the profile analysis method to further remove large background objects without affecting bridges. According to profile analysis, points on the same object usually have similar heights, whereas differing objects, especially tall objects and the ground, have major height differences. Thus, by appropriately setting the neighborhood threshold and height threshold, points belonging to the same object can be easily connected and removed. In particular, as both ends of a bridge are always smoothly connected to the ground, the bridges belong to the same connected domain with ground. Consequently, potential bridges will remain when other large objects are removed. As shown in Fig. 4, the process of large objects removal based on profile analysis can be summarized in the following steps:
Fig. 6 shows a simple case of profile analysis. In the figure, the Lidar data are partitioned into three profiles in the directions 45°, 90°, 135°. The final connected domains are derived by overlaying and connecting stripes in the three profiles. Owing to the significant amount of data in Lidar, generating the MST in step (2) and ascertaining all the connected domains in step (4) is time-consuming. Conventional methods, such as the Prim method, are
(1) Partition the Lidar data into four profiles, which are used to form stripes in four directions. Usually the four directions are chosen as 0°, 45°, 90°, 135° according to coordinate rotation theory, and the stride width in each profile is set to be slightly larger than the average point space. (2) For each profile in a specified direction, stripes are formed by
Fig. 3. Vegetation removal. (a) Lidar data. (b) Vegetation removal result after using the proposed 3D opening operator. 3
Optics and Laser Technology 121 (2020) 105790
L. Li, et al.
Fig. 4. Large objects removal process.
where h is the height of the current point and hi is the height of the ith neighborhood point, for i = 1, 2, … , n. Usually the number of neighborhood points n is set as 10. Further, the smoothness of a point is defined as the standard deviation of the differences between the actual surface and the surface fitted by its m neighborhood points with the formula
often inefficient and result in insufficient-memory problems. To avoid these problems and optimize the two processes, the proposed method utilizes union-find sets to quickly merge sets and decide the set to which an element belongs. Conventionally, union-find sets is realized by first initializing each element as the ancestor node of a new set, then deciding which set each element belongs to and determining its ancestor node. Subsequently, the two sets are merged by assigning the ancestor of one set to the other’s ancestor. To further optimize this process, the path compression method is used. As shown in Fig. 7, once the ancestor node is found by recursion, through pointing all of its children modes to itself during backtracking, the complexity of element search in the following steps decreases from O(n) to O(1). In addition, merging sets based on rank, specifically, merging the set with fewer elements into the one with more elements, reduces the height of the tree.
xsmoothness =
i=1
(5)
1. Calculate the two quantitative measurements: discontinuity and smoothness. Because points on the long sides of the bridge that have large height differences to the ground and points on the short sides that are smoothly connected to the ground must exist, their discontinuity and smoothness are used to select the points belonging to bridges. Usually, the discontinuity threshold is set as the height of the bridge, such as 3 m, while the smoothness threshold is set as 0.25. 2. Partition the Lidar points into profiles in six directions (chosen as 0°, 30°, 60°, 90°, 120°, 150° according to coordinate rotation theory), as in the large objects removal process. 3. For each profile in a specified direction, according to the given slope (such as 45°) and height threshold (such as three times the average point space), connect the points and assign the shapes of the lines (see Fig. 5) to their end points. 4. Overlap the six profiles, and select the seed points. Because a point in different profiles may have different shapes, the seed points are
This section describe how the bridges are extracted from the remaining Lidar data based on the profile analysis method discussed above and the characteristics of bridges. It is well-known that the both ends of a bridge are connected to the ground and the two sides have a large height difference with ground. Regarding a bridge as a rectangle, there are usually many points in its two longer sides whose height is much greater than 3 m (for practical purposes, the bridge height is more than 3 m), and the other two sides are smoothly connected to the ground and are longer than 3 m (that is, the bridge width is more than 3 m). In order to quantitatively measure the discontinuity and smoothness, two parameters, x discontinuity and x smoothness are defined. The discontinuity of a point is the maximum height difference between the current point and its neighborhood points, as shown in Eq. (4), i = 1,2,... n
m
∑ (Δhi − μ)2
where hi represents the height of the ith neighborhood point, and its corresponding point in the fitted plane is i’ with height hi’, △hi = hi hi′, µ is the mean value of △hi. Usually m is set as 20, and the fitted surface is obtained based on the least-squares fitting method. With the above characteristics of bridges, as shown in Fig. 8, bridge extraction based on profile analysis is realized using the following steps:
2.4. Bridge extraction
x discontinuity = max {h − hi}
1 m
(4)
Fig. 5. Shapes of lines and points. (Set domain ① as the reference standard.) (a) no shape; (b) lowered left; (c) raised left; (d) terraced; (e) lowered; (f) raised. 4
Optics and Laser Technology 121 (2020) 105790
L. Li, et al.
Fig. 6. Profile analysis with regular points for expression.
Fig. 7. Optimization of Union-find Sets using path compression.
Fig. 8. Bridge detection.
3. Parallel implementation with OpenMP
selected as those that have the most significant shape raise. 5. Construct Delaunay triangulation using seed points, and delete those edges that do not meet the neighborhood and slope threshold (for example, 2.5 times average point space and 30°). Determine all the connected domains that could be potential bridges.
Because data processing based on profile analysis in different profiles is irrelevant, parallelization can be applied to further improve detection efficiency for every profile. As a library that supports parallel design with shared memory, OpenMP is especially suitable for multicore CPU parallel programming [19]. It uses the standard parallel mode fork-join of shared storage, as shown in Fig. 9. When the program begins, only the main thread exists and operations are sequential. Subsequently, at the derived point, several parallel threads used for parallel data processing are produced. After the parallel threads end, the main thread starts again at the merge
As the size of the connected domain for a bridge is limited, which means that the number of points representing a bridge is limited with the known average point space, set an interval number to remove some areas that are obviously not bridges. Finally, accurately extract the bridges with the characteristics of bridges described above.
5
Optics and Laser Technology 121 (2020) 105790
L. Li, et al.
Fig. 9. Fork-Join. Table 1 Comparison of computation time. Data size (number of Lidar points)
22259 74317 173598 486800
Sithole and Vosselman’s method
Our proposed method
Time used in large objects removal process (s)
Time used in bridge extraction process (s)
Time used in large objects removal process (s)
Time used in bridge extraction process (s)
31 103 196 435
5 12 40 65
10 30 60 124
2 3 7 16
Fig. 10. Bridge detection result for the first group of discrete point cloud data. (a) Original point clouds data; (b) DEM; (c) vegetation removal; (d) large objects removal; (e) bridges extraction; (f) bridges shown in original data.
different sizes at the same position of the dataset. We selected four directions (0◦, 45◦, 90◦, 135◦) int the large objects removal process and another six directions (0◦, 30◦, 60◦, 90◦, 120◦, 150◦) in the bridge extraction process. Compared to the computation time required by Sithole and Vosselman’s method, it is clear that by applying union-find sets and OpenMP our method reduces the computation time by almost 60%. To intuitively show the accuracy of the proposed method for bridge detection, three groups of representative discrete point cloud data were chosen. The first group was urban point cloud data with average point space of 0.75 m and 98,696 points, containing vegetation, roads, bridges, buildings, and other objects. The second one contained vegetation, roads, and bridges, with 95,131 points and average point space
point. Based on these features of OpenMP, each profile in the large objects removal and bridge extraction parts can be processed by one parallel thread, that is, there are four and six parallel threads in these two parts, respectively.
4. Results To demonstrate the efficiency and accuracy of the proposed method, experimental data were selected from a complex point dataset that contains vegetation, buildings, bridges, etc. There are 486,800 points in this dataset, and the average point space is 0.8 m. Table 1 shows the experiment results based on four sub-graphs derived by cutting out 6
Optics and Laser Technology 121 (2020) 105790
L. Li, et al.
Fig. 11. Bridge detection result for the second group of discrete point cloud data. (a) Original point clouds data; (b) DEM; (c) vegetation removal; (d) large objects removal; (e) bridges extraction; (f) bridges shown in original data.
Fig. 12. Bridge detection result for the third group of discrete point cloud data. (a) Original point clouds data; (b) DEM; (c) vegetation removal; (d) large objects removal; (e) bridges extraction; (f) bridges shown in original data.
Figs. 10(f), 11(f) and 12(f) show the bridge extraction result in the original point cloud data.
of 0.8 m. The second group contained more vegetation and the bridge location was less obvious than the first one. The third group was mountain area point cloud data points and average point space of 1.59 m and 11,837,440 points, containing vegetation, roads, bridges, mountain and other objects. Figs. 10, 11 and 12 present the experimental results corresponding to the above three groups of discrete point cloud data. In each figure, (a) is the original point cloud data and (b) shows the DEM description of the point cloud data, from which it can be seen that the data contain many buildings and vegetation, as well as a small number of small objects and bridges. Figs. 10(c), 11(c) and 12(c) are the result of the last and single echo filtering, together with the process of the morphological opening operation. Subsequently, the vegetation and some small objects can be effectively removed. Figs. 10(d), 11(d) and 12(d) show the large objects removal results, the buildings and other large objects have been removed and the potential bridges have been retained. Figs. 10(e), 11(e) and 12(e) give the extraction result of the bridge point.
5. Conclusions This paper proposed an improved detection method based on morphology and profile analysis for bridge extraction from Lidar data. In the proposed method, preprocessing by denoising and selecting the last and single echo data can filter a certain amount of vegetation and noise points. The morphological operations based on a newly designed morphological operator for 3D discrete points can not only effectively remove most of the vegetation and some small objects, but also avoid the detrimental effects caused by vegetation. Considering that direct calculation of the neighborhood of the current point needed in morphological operations is time-consuming, a gridding method that significantly reduces the calculation time is introduced. Using profile analysis, the large objects in the background are removed, and by 7
Optics and Laser Technology 121 (2020) 105790
L. Li, et al.
combining profile analysis and the topological characteristics of bridges further, the bridges are extracted successfully. The profile analysis method used is optimized by applying union-find sets to speed up the process of generating the MST and determining the connected domains. Moreover, further improvement in computation efficiency is achieved by parallel implementation with OpenMP. The proposed method is independent of the bridge shape also. Experimental results show that bridges can be extracted from Lidar data with high efficiency and accuracy. But for the case of bridges across water, because little or no echo return to the light collector owing to mirror reflection, the proposed method cannot be able to accurately extract them. Further work will focus on solving this problem.
Model, Remote Sens. Environ. 121 (2012) 224–235. [5] F. Rottensteiner, J. Trinder, S. Clode, K. Kubik, Using the Dempster-Shafer Method for the fusion of Lidar data and multi-spectral images for building detection, Inform. Fusion 6 (2005) 283–300. [6] U.C. Benz, P. Hofmann, G. Willhauck, I. Lingenfelder, M. Heynen, Multi-resolution, object-oriented fuzzy analysis of remote sensing data for gis-ready information, ISPRS J. Photogram. Remote Sens. 58 (2004) 239–258. [7] B.N. Jyothi, G.R. Babu, I.V.M. Krishna, Object oriented and multi-scale image analysis: strengths, weaknesses, opportunities and threats-A, review, J. Comput. Sci. 4 (2008) 706–712. [8] M. Baatz, Multiresolution segmentation-an optimization approach for high quality multi-scale image segmentation, Beiträge zum AGIT-Symp. (2000). [9] S.S. Hwang, J.K. Aggarwal, Detection and segmentation of man-made objects in outdoor scenes: concrete bridges, J. Opt. Soc. Am. A 6 (1989) 938–950. [10] An Automatic Bridge Detection Technique for Multispectral Images, IEEE Trans. Geosci. Remote Sens. 46 (2008) 2720–2727. [11] G. Sithole, G. Vosselman, Filtering of Airborne Laser Scanner Data Based on Segmented Point Clouds, 2005. [12] G. Sithole, G. Vosselman, Bridge detection in airborne laser scanner data, Isprs J. Photogram. Remote Sens. (2007) 61–71. [13] L.I. Yunfan, Extracting bridges from airborne lidar data based on terrain features, Geomat. Inform. Sci. Wuhan Univ. 36 (2011) 552–555. [14] A. Samberg, An implementation of the Asprs Las Standard, ISPRS – Int. Arch. Photogram., Remote Sens. Spatial Inform. Sci., 2007, pp. 3–52. [15] Filtering Airborne Laser Scanning Data with Morphological Methods, Photogram. Eng. Remote Sens. 73 (2015) 175–185. [16] A Gradient-Constrained Morphological Filtering Algorithm for Airborne Lidar, Opt. Laser Technol. 54 (2013) 288–296. [17] A.K. Singh, Improved hybrid algorithm for robust and imperceptible multiple watermarking using digital images, Multimed. Tools Appl. 76 (6) (2017) 8881–8900, https://doi.org/10.1007/s11042-016-3514-z. [18] G. Sithole, Filtering of Laser Altimetry Data Using Slope Adaptive Filter Vol. 22 Annapolis, Maryland, 2001. [19] Serial and Parallel Implementation of Shortest Path Algorithm in the Optimization of Public Transport Travel, 2011.
Acknowledgements This study is funded by the National Key Research and Development Program of China (2016YFB0501403), the Scientific and Technological Leading Talent Fund of the National Administration of Surveying, Mapping and Geo-information (2014). References [1] M. Awrangjeb, C. Zhang, F.C.S, Automatic extraction of building roofs using lidar data and multispectral imagery, Isprs J. Photogram. Remote Sens. 83 (2013) 1–18. [2] G. Sohn, I. Downman, Data fusion of high-resolution satellite imagery and Lidar data for automatic building extraction A, Isprs J. Photogram. Remote Sens. 62 (2007) 43–63. [3] P. Axelsson, Dem generation from laser scanner data using adaptive tin models, Int. Arch. Photogram. Remote Sens. (2000) 33. [4] Accuracy Assessment and Correction of a Lidar-Derived Salt Marsh Digital Elevation
8
Contents lists available at ScienceDirect
Optics and Laser Technology journal homepage: www.elsevier.com/locate/optlastec
Full length article
An improved detection method based on morphology and profile analysis for bridge extraction from Lidar
T
⁎
Lin Lia, Wei Rongb,a, , Fei Sub, Xiaoyu Xingb a b
Collaborative Innovation Centre of Geospatial Technology, Wuhan University, 129 Luoyu Road, Wuhan 430079, China School of Resource and Environmental Sciences, Wuhan University, 129 Luoyu Road, Wuhan 430079, China
H I GH L IG H T S
dimensional discrete points morphological operator, removing vegetation. • Three sets incorporated with the profile analysis method, removing large background objects. • Union-find optimized profile analysis method with the topological characteristics, extracting the bridges. • The • OpenMP to improve the computation efficiency.
A R T I C LE I N FO
A B S T R A C T
Keywords: Bridge extraction Three dimensional morphology Profile analysis OpenMP
Extraction of bridges from light detection and ranging (Lidar) images is a difficult problem with low detection accuracy and detection efficiency because of strong dependence on bridge shapes, the influence of vegetation and the large amount of data. This paper proposes an improved method based on morphology and profile analysis to extract bridges with removing background objects such as vegetation. To remove vegetation effectively, a new morphological three dimensional (3D) discrete points operator is applied which the gridding method is utilized to obtain the neighborhoods of the discrete points. Subsequently, union-find sets is incorporated with the profile analysis method to optimize the process of generating the minimum spanning tree (MST) and determining connected domains, which are used to remove large objects and leave the bridges behind. By combining the optimized profile analysis method with the topological characteristics of bridges, bridges are extracted without dependence on their shapes. Finally, in order to improve the computation efficiency of the Lidar data, the OpenMP is employed. Experimental results show that the proposed method can extract the bridges from Lidar data effectively.
1. Introduction Light detection and ranging (Lidar) is a new airborne remote sensing technology that can be used to quickly and accurately obtain the position and elevation, i.e., the point cloud data, of ground objects. Currently, the main areas in which Lidar is applied include aspects of building extraction and reconstruction [1,2], digital elevation model (DEM) generation [3,4], and residential area recognition based on Lidar and remote sensing images [5]. The bridge not only play an important role in city planning but are also is a crucial military targets in war. Consequently, accurate and rapid detection of bridge has important research value. Most existing bridge detection methods are based on optical remote sensing aerospace images instead of Lidar data. These methods can be classified into
⁎
object-oriented analysis methods [6–8], and knowledge-based methods [9,10]. These methods were developed for the case that the bridges are over water, however, the bridges cannot efficiently extract bridges on land in general. This paper proposes an improved method based on Sithole and Vosselman’s algorithm that overcomes the problems outlined above. The proposed method gradually removes noise points, vegetation, and large objects to extract the bridges from Lidar images accurately. The gridding method is applied for removing the vegetation quickly by a newly designed 3D morphological operator. The profile analysis method is optimized using union-find sets and topology characteristics, through the optimized method, the large background objects can be removed and the bridges are extracted accurately. In order to improve the efficiency of the method further, the parallel processing with
Corresponding author. E-mail addresses: [email protected] (W. Rong), [email protected] (F. Su), [email protected] (X. Xing).
https://doi.org/10.1016/j.optlastec.2019.105790 Received 5 November 2018; Received in revised form 8 May 2019; Accepted 25 August 2019 0030-3992/ © 2019 Elsevier Ltd. All rights reserved.
Optics and Laser Technology 121 (2020) 105790
L. Li, et al.
Fig. 1. Process flow of the proposed method of bridge extraction in this paper.
each other. As a common image processing theory, mathematical morphology is usually used for ground point extraction and point cloud filtering in Lidar data processing [15,16]. In this part of the process, the proposed 3D morphological opening operator can be utilized to remove most of the vegetation and some small objects without affecting large objects. Assuming that A represents a point set from 2D images, for each point a, where a ∈ A, the structure element B is defined as the set of neighborhood points of a. The morphological opening operation is carried out by first using B to corrode A and then expanding the corrosion result with B, as shown in Eq. (1) [17]:
OpenMP is used to obtain the stripe information from various directions in profile analysis. 2. Methods The proposed bridge extraction method consists of four main steps: preprocessing, vegetation removal, large objects removal, and bridge detection as shown in Fig. 1. Following preprocessing with denoising and filtering, a new 3D discrete points morphological operator is applied to remove vegetation. Next, large background objects are removed based on the optimized profile analysis method. Finally, combined with the topology characteristics of bridges, the optimized profile analysis method is used to extract the bridges.
(1)
A°B = (AΘB ) ⊕ B
To apply opening operations on Lidar data, we redefine A as the set of point cloud data of Lidar. For each point a, where a ∈ A, set B is the neighborhood points of a in 3D data space, and a.z represents the height of point a. The corrosion operation on A assigns the height of the lowest point in B to a.z, while the expansion operation on A assigns the height of the highest point in B to a.z. The opening operation on A can be described as:
2.1. Preprocessing Lidar data are obtained from the Lidar system, which comprises a global positioning system (GPS), scanner, inertial measurement unit (IMU), and imaging devices. In general, it uses a laser light to scan the earth and other objects, by calculating the echo time, the equipment can obtain detailed points information about the surface of the ground or the objects, such as height and coordinates. However, the laser light may strike special objects, such as dust in the air, glasses, and metallic shimmering surfaces. In such situations, reflection and refraction of the light will happen. The points of these special data are defined as noise points. Further, owing to the strong penetrability of laser light, the redundant information caused by multiple reflections may exist in Lidar data also. So, it is necessary to filter the noise points and redundant information from the original Lidar data before subsequent operations.
AΘB = \{ a ∈ A|a. z = min (B. z )\} A ⊕ B = \{ a ∈ A|a. z = max (B. z )\} A°B = (AΘB ) ⊕ B
(2)
In general, the neighborhood of one point can be easily obtained in 2D images according to the coordinates. In contrast, in Lidar data, owing to the huge amount of 3D discrete points, the tremendous calculations on the distances between the current point and all the other points significantly reduces the speed of operations. Thus, the gridding method is utilized to rapidly obtain the neighborhood of 3D discrete points. Instead of directly calculating the neighborhood of the current point, by first gridding the Lidar data and ascertaining the neighborhood of the grid containing the current point, the neighborhood of the current point can be determined from the current grid and its neighborhood. The procedure is as follows:
2.1.1. Noise filtering The noise points derived from dangling objects often cause some abnormal data, in which the height differences between the data collector and the ground or the objects are greater than the real value. A simple and efficient way to remove these noise points is to set an appropriate height threshold. For example, if the real height difference is 500 m, we can assign 500 as the threshold, which means data greater than 500 will be discarded.
1. Grid the discrete points (Fig. 2(b)) with grid width twice the average point space, to ensure that all the neighborhood points (usually eight points) are included in the current grid and its neighborhood; here the current grid is defined as the grid containing the current point (point 0). 2. Calculate the distances between all the other points in the nine grids and the current point, and select the eight nearest points using the Top-K sorting algorithm.
2.1.2. Single and last echo data filtering As the penetrability of laser light is very high, when the light reaches an object, especially vegetation, multiple reflections occur. In these reflections, the last echo data usually come from ground. So, selecting the last and single echo data can filtering a certain amount of vegetation and noise points, and also reduce the volume of Lidar data. However, buildings and other large objects do not produce multiple reflections because of their special materials, so this step has no effect on the points representing them. Based on the Lidar standard data format (LAS) [14], single and last echo data filtering can be used to extract those points whose echo numbers are equal to their echo times. After preprocessing, the Lidar data only consist of the points representing the remaining vegetation, large objects, and ground.
This method significantly reduces the search scope and the amount of calculation. Note that the average point space mentioned in step (1) is calculated based on the data format of Lidar points using Eq. (3):
AveragePointSpacing =
2.2. Vegetation removal
(max X − min X )(max Y − min Y ) PointsNum
(3)
where maxX, minX, maxY, and minY are parameters of the data scope, and PointsNum is the overall number of points. Fig. 3 gives an example of vegetation removal with the above 3D morphological opening operations. In the figure, the Lidar data include vegetation, buildings, and roads. The number of points is 98696, and the average point space is 0.75 m. The result in Fig. 3(b) shows that the proposed 3D opening operator removed a reasonable amount of vegetation and some other small objects.
In Lidar data, the points on the same large object have virtually the same height, except for the boundary of the object—especially the area connected with ground, where the points usually have large height differences from the ground. In other words, the points on one large object exist as “patches.” As the ground is usually next to vegetation or some small objects, the points often have greater distances between 2
Optics and Laser Technology 121 (2020) 105790
L. Li, et al.
Fig. 2. Procedure for obtaining the neighborhood of the discrete points using the gridding method: (a) Discrete Lidar Points; (b) Points Gridded; (c) Eight Neighborhood Points (1–8) of Point 0.
constructing Delaunay triangulations and generating the minimum spanning tree (MST). The height threshold and neighborhood threshold, generally set as one-third of average point space and twice the average point space, respectively, are used to delete the edges with high weights in the MST. (3) Overlap the four profiles, which results in some stripes in different directions being connected into a larger connected domain. (4) Determine all the connected domains, each of which represents an object. To accurately indicate the feature of lines and points in the connected domain, six different shapes of lines and points are defined [12], as shown in Fig. 5. For example, if connected region ① only has one adjacent connected domain on the left or right side, and ① is the lower region, then for the shapes of connected region ①, all the lines and points in ① are defined as lowered. Based on the shape of the lines and points in each connected domain [11], set the threshold to determine whether the connected domain belongs to a large object. If it does, then all the points in that connected domain will be directly removed. (5) Iterate steps (1)–(4) between one and four times, depending on the actual components of the point cloud data.
2.3. Large objects removal After filtering the noise points and small objects, the remaining Lidar data mainly consist of the points representing the large objects and the ground. Although existing methods based on the geometric characteristics of Lidar points [3,13,18] perform well in extracting ground points to generate a DEM, they often remove bridges as well when removing large objects. In contrast, the proposed method utilizes the profile analysis method to further remove large background objects without affecting bridges. According to profile analysis, points on the same object usually have similar heights, whereas differing objects, especially tall objects and the ground, have major height differences. Thus, by appropriately setting the neighborhood threshold and height threshold, points belonging to the same object can be easily connected and removed. In particular, as both ends of a bridge are always smoothly connected to the ground, the bridges belong to the same connected domain with ground. Consequently, potential bridges will remain when other large objects are removed. As shown in Fig. 4, the process of large objects removal based on profile analysis can be summarized in the following steps:
Fig. 6 shows a simple case of profile analysis. In the figure, the Lidar data are partitioned into three profiles in the directions 45°, 90°, 135°. The final connected domains are derived by overlaying and connecting stripes in the three profiles. Owing to the significant amount of data in Lidar, generating the MST in step (2) and ascertaining all the connected domains in step (4) is time-consuming. Conventional methods, such as the Prim method, are
(1) Partition the Lidar data into four profiles, which are used to form stripes in four directions. Usually the four directions are chosen as 0°, 45°, 90°, 135° according to coordinate rotation theory, and the stride width in each profile is set to be slightly larger than the average point space. (2) For each profile in a specified direction, stripes are formed by
Fig. 3. Vegetation removal. (a) Lidar data. (b) Vegetation removal result after using the proposed 3D opening operator. 3
Optics and Laser Technology 121 (2020) 105790
L. Li, et al.
Fig. 4. Large objects removal process.
where h is the height of the current point and hi is the height of the ith neighborhood point, for i = 1, 2, … , n. Usually the number of neighborhood points n is set as 10. Further, the smoothness of a point is defined as the standard deviation of the differences between the actual surface and the surface fitted by its m neighborhood points with the formula
often inefficient and result in insufficient-memory problems. To avoid these problems and optimize the two processes, the proposed method utilizes union-find sets to quickly merge sets and decide the set to which an element belongs. Conventionally, union-find sets is realized by first initializing each element as the ancestor node of a new set, then deciding which set each element belongs to and determining its ancestor node. Subsequently, the two sets are merged by assigning the ancestor of one set to the other’s ancestor. To further optimize this process, the path compression method is used. As shown in Fig. 7, once the ancestor node is found by recursion, through pointing all of its children modes to itself during backtracking, the complexity of element search in the following steps decreases from O(n) to O(1). In addition, merging sets based on rank, specifically, merging the set with fewer elements into the one with more elements, reduces the height of the tree.
xsmoothness =
i=1
(5)
1. Calculate the two quantitative measurements: discontinuity and smoothness. Because points on the long sides of the bridge that have large height differences to the ground and points on the short sides that are smoothly connected to the ground must exist, their discontinuity and smoothness are used to select the points belonging to bridges. Usually, the discontinuity threshold is set as the height of the bridge, such as 3 m, while the smoothness threshold is set as 0.25. 2. Partition the Lidar points into profiles in six directions (chosen as 0°, 30°, 60°, 90°, 120°, 150° according to coordinate rotation theory), as in the large objects removal process. 3. For each profile in a specified direction, according to the given slope (such as 45°) and height threshold (such as three times the average point space), connect the points and assign the shapes of the lines (see Fig. 5) to their end points. 4. Overlap the six profiles, and select the seed points. Because a point in different profiles may have different shapes, the seed points are
This section describe how the bridges are extracted from the remaining Lidar data based on the profile analysis method discussed above and the characteristics of bridges. It is well-known that the both ends of a bridge are connected to the ground and the two sides have a large height difference with ground. Regarding a bridge as a rectangle, there are usually many points in its two longer sides whose height is much greater than 3 m (for practical purposes, the bridge height is more than 3 m), and the other two sides are smoothly connected to the ground and are longer than 3 m (that is, the bridge width is more than 3 m). In order to quantitatively measure the discontinuity and smoothness, two parameters, x discontinuity and x smoothness are defined. The discontinuity of a point is the maximum height difference between the current point and its neighborhood points, as shown in Eq. (4), i = 1,2,... n
m
∑ (Δhi − μ)2
where hi represents the height of the ith neighborhood point, and its corresponding point in the fitted plane is i’ with height hi’, △hi = hi hi′, µ is the mean value of △hi. Usually m is set as 20, and the fitted surface is obtained based on the least-squares fitting method. With the above characteristics of bridges, as shown in Fig. 8, bridge extraction based on profile analysis is realized using the following steps:
2.4. Bridge extraction
x discontinuity = max {h − hi}
1 m
(4)
Fig. 5. Shapes of lines and points. (Set domain ① as the reference standard.) (a) no shape; (b) lowered left; (c) raised left; (d) terraced; (e) lowered; (f) raised. 4
Optics and Laser Technology 121 (2020) 105790
L. Li, et al.
Fig. 6. Profile analysis with regular points for expression.
Fig. 7. Optimization of Union-find Sets using path compression.
Fig. 8. Bridge detection.
3. Parallel implementation with OpenMP
selected as those that have the most significant shape raise. 5. Construct Delaunay triangulation using seed points, and delete those edges that do not meet the neighborhood and slope threshold (for example, 2.5 times average point space and 30°). Determine all the connected domains that could be potential bridges.
Because data processing based on profile analysis in different profiles is irrelevant, parallelization can be applied to further improve detection efficiency for every profile. As a library that supports parallel design with shared memory, OpenMP is especially suitable for multicore CPU parallel programming [19]. It uses the standard parallel mode fork-join of shared storage, as shown in Fig. 9. When the program begins, only the main thread exists and operations are sequential. Subsequently, at the derived point, several parallel threads used for parallel data processing are produced. After the parallel threads end, the main thread starts again at the merge
As the size of the connected domain for a bridge is limited, which means that the number of points representing a bridge is limited with the known average point space, set an interval number to remove some areas that are obviously not bridges. Finally, accurately extract the bridges with the characteristics of bridges described above.
5
Optics and Laser Technology 121 (2020) 105790
L. Li, et al.
Fig. 9. Fork-Join. Table 1 Comparison of computation time. Data size (number of Lidar points)
22259 74317 173598 486800
Sithole and Vosselman’s method
Our proposed method
Time used in large objects removal process (s)
Time used in bridge extraction process (s)
Time used in large objects removal process (s)
Time used in bridge extraction process (s)
31 103 196 435
5 12 40 65
10 30 60 124
2 3 7 16
Fig. 10. Bridge detection result for the first group of discrete point cloud data. (a) Original point clouds data; (b) DEM; (c) vegetation removal; (d) large objects removal; (e) bridges extraction; (f) bridges shown in original data.
different sizes at the same position of the dataset. We selected four directions (0◦, 45◦, 90◦, 135◦) int the large objects removal process and another six directions (0◦, 30◦, 60◦, 90◦, 120◦, 150◦) in the bridge extraction process. Compared to the computation time required by Sithole and Vosselman’s method, it is clear that by applying union-find sets and OpenMP our method reduces the computation time by almost 60%. To intuitively show the accuracy of the proposed method for bridge detection, three groups of representative discrete point cloud data were chosen. The first group was urban point cloud data with average point space of 0.75 m and 98,696 points, containing vegetation, roads, bridges, buildings, and other objects. The second one contained vegetation, roads, and bridges, with 95,131 points and average point space
point. Based on these features of OpenMP, each profile in the large objects removal and bridge extraction parts can be processed by one parallel thread, that is, there are four and six parallel threads in these two parts, respectively.
4. Results To demonstrate the efficiency and accuracy of the proposed method, experimental data were selected from a complex point dataset that contains vegetation, buildings, bridges, etc. There are 486,800 points in this dataset, and the average point space is 0.8 m. Table 1 shows the experiment results based on four sub-graphs derived by cutting out 6
Optics and Laser Technology 121 (2020) 105790
L. Li, et al.
Fig. 11. Bridge detection result for the second group of discrete point cloud data. (a) Original point clouds data; (b) DEM; (c) vegetation removal; (d) large objects removal; (e) bridges extraction; (f) bridges shown in original data.
Fig. 12. Bridge detection result for the third group of discrete point cloud data. (a) Original point clouds data; (b) DEM; (c) vegetation removal; (d) large objects removal; (e) bridges extraction; (f) bridges shown in original data.
Figs. 10(f), 11(f) and 12(f) show the bridge extraction result in the original point cloud data.
of 0.8 m. The second group contained more vegetation and the bridge location was less obvious than the first one. The third group was mountain area point cloud data points and average point space of 1.59 m and 11,837,440 points, containing vegetation, roads, bridges, mountain and other objects. Figs. 10, 11 and 12 present the experimental results corresponding to the above three groups of discrete point cloud data. In each figure, (a) is the original point cloud data and (b) shows the DEM description of the point cloud data, from which it can be seen that the data contain many buildings and vegetation, as well as a small number of small objects and bridges. Figs. 10(c), 11(c) and 12(c) are the result of the last and single echo filtering, together with the process of the morphological opening operation. Subsequently, the vegetation and some small objects can be effectively removed. Figs. 10(d), 11(d) and 12(d) show the large objects removal results, the buildings and other large objects have been removed and the potential bridges have been retained. Figs. 10(e), 11(e) and 12(e) give the extraction result of the bridge point.
5. Conclusions This paper proposed an improved detection method based on morphology and profile analysis for bridge extraction from Lidar data. In the proposed method, preprocessing by denoising and selecting the last and single echo data can filter a certain amount of vegetation and noise points. The morphological operations based on a newly designed morphological operator for 3D discrete points can not only effectively remove most of the vegetation and some small objects, but also avoid the detrimental effects caused by vegetation. Considering that direct calculation of the neighborhood of the current point needed in morphological operations is time-consuming, a gridding method that significantly reduces the calculation time is introduced. Using profile analysis, the large objects in the background are removed, and by 7
Optics and Laser Technology 121 (2020) 105790
L. Li, et al.
combining profile analysis and the topological characteristics of bridges further, the bridges are extracted successfully. The profile analysis method used is optimized by applying union-find sets to speed up the process of generating the MST and determining the connected domains. Moreover, further improvement in computation efficiency is achieved by parallel implementation with OpenMP. The proposed method is independent of the bridge shape also. Experimental results show that bridges can be extracted from Lidar data with high efficiency and accuracy. But for the case of bridges across water, because little or no echo return to the light collector owing to mirror reflection, the proposed method cannot be able to accurately extract them. Further work will focus on solving this problem.
Model, Remote Sens. Environ. 121 (2012) 224–235. [5] F. Rottensteiner, J. Trinder, S. Clode, K. Kubik, Using the Dempster-Shafer Method for the fusion of Lidar data and multi-spectral images for building detection, Inform. Fusion 6 (2005) 283–300. [6] U.C. Benz, P. Hofmann, G. Willhauck, I. Lingenfelder, M. Heynen, Multi-resolution, object-oriented fuzzy analysis of remote sensing data for gis-ready information, ISPRS J. Photogram. Remote Sens. 58 (2004) 239–258. [7] B.N. Jyothi, G.R. Babu, I.V.M. Krishna, Object oriented and multi-scale image analysis: strengths, weaknesses, opportunities and threats-A, review, J. Comput. Sci. 4 (2008) 706–712. [8] M. Baatz, Multiresolution segmentation-an optimization approach for high quality multi-scale image segmentation, Beiträge zum AGIT-Symp. (2000). [9] S.S. Hwang, J.K. Aggarwal, Detection and segmentation of man-made objects in outdoor scenes: concrete bridges, J. Opt. Soc. Am. A 6 (1989) 938–950. [10] An Automatic Bridge Detection Technique for Multispectral Images, IEEE Trans. Geosci. Remote Sens. 46 (2008) 2720–2727. [11] G. Sithole, G. Vosselman, Filtering of Airborne Laser Scanner Data Based on Segmented Point Clouds, 2005. [12] G. Sithole, G. Vosselman, Bridge detection in airborne laser scanner data, Isprs J. Photogram. Remote Sens. (2007) 61–71. [13] L.I. Yunfan, Extracting bridges from airborne lidar data based on terrain features, Geomat. Inform. Sci. Wuhan Univ. 36 (2011) 552–555. [14] A. Samberg, An implementation of the Asprs Las Standard, ISPRS – Int. Arch. Photogram., Remote Sens. Spatial Inform. Sci., 2007, pp. 3–52. [15] Filtering Airborne Laser Scanning Data with Morphological Methods, Photogram. Eng. Remote Sens. 73 (2015) 175–185. [16] A Gradient-Constrained Morphological Filtering Algorithm for Airborne Lidar, Opt. Laser Technol. 54 (2013) 288–296. [17] A.K. Singh, Improved hybrid algorithm for robust and imperceptible multiple watermarking using digital images, Multimed. Tools Appl. 76 (6) (2017) 8881–8900, https://doi.org/10.1007/s11042-016-3514-z. [18] G. Sithole, Filtering of Laser Altimetry Data Using Slope Adaptive Filter Vol. 22 Annapolis, Maryland, 2001. [19] Serial and Parallel Implementation of Shortest Path Algorithm in the Optimization of Public Transport Travel, 2011.
Acknowledgements This study is funded by the National Key Research and Development Program of China (2016YFB0501403), the Scientific and Technological Leading Talent Fund of the National Administration of Surveying, Mapping and Geo-information (2014). References [1] M. Awrangjeb, C. Zhang, F.C.S, Automatic extraction of building roofs using lidar data and multispectral imagery, Isprs J. Photogram. Remote Sens. 83 (2013) 1–18. [2] G. Sohn, I. Downman, Data fusion of high-resolution satellite imagery and Lidar data for automatic building extraction A, Isprs J. Photogram. Remote Sens. 62 (2007) 43–63. [3] P. Axelsson, Dem generation from laser scanner data using adaptive tin models, Int. Arch. Photogram. Remote Sens. (2000) 33. [4] Accuracy Assessment and Correction of a Lidar-Derived Salt Marsh Digital Elevation
8