## 2 Data and Methods
### LiDAR Feature Extraction  
LiDAR is an active sensing system that operates by measuring the round-trip timing of a laser beam from an object in order to determine the distance from the sensor to the target. By analyzing the laser time range in combination with the scan angle and spatial coordinates of the laser scanner, the spatial coordinates of an object as (X, Y, Z) are obtained. Along with spatial data, intensity values, number of returns, point classification values, and GPS times are recorded. However, the complete geometric properties of the targetted objects in the point cloud cannot be solely represented by their spatial information (X, Y, and Z), though they can be characterized by their geometric features. In order to avoid the heavy computational load of using complex point cloud data, semantic information, extracted by generating geometric features, can be used instead. As presented in Figure 2, feature extraction was carried out in four steps. In the first step, points whose noise levels are high were eliminated from the LiDAR point cloud. After noise removal, outlier detection was performed by removing the points whose average distance from their neighbouring points exceeded a given threshold value. In the next step, the statistical relationship of the points was formed by creating a new data structure based on neighbouring points, similar to the second step. But this time, the neighbouring points were determined using a clean point cloud by applying a k-dimensional tree (k-d tree) algorithm in 3D. In the final step, eigenvalues and the 3D geometric features were calculated. There are various methods available for determining the neighborhood of points in LiDAR data. The neighborhood of points can be constructed using a spherical radius or parameterized based on the number of closest neighbours in 2D or 3D space (i.e., k-NN methods). In this study, k-NN based neighbourhood selection was constructed via the well-known \(k\)-d tree algorithm, in which the data are partitioned using a binary tree  
Figure 1: Satellite image (**a**), rendering of point cloud (**b**) and their overlapping image composite (**c**), and comparison of satellite image and LiDAR data side by side for the area of interest enclosed by the black box (**d**). Figure 2: Feature extraction pipeline. structure [37]. As a result of the neighbourhood selection, each query point (\(P(i)\)) and its neighbors were indexed. The distance between the query point and its \(k\) neighbours was determined and indexed in order from nearest to farthest. In the first part of the feature extraction, the features that represent the geometric 3D properties of the neighbourhood were calculated using the height of each query point, i.e., the absolute height \(H\), the height difference between the query point and its neighbouring points \(\Delta H_{\text{kNN,3D}}\), the standard deviation of absolute height within neighbouring points \(\sigma H_{\text{kNN,3D}}\), local 3D neighbourhood represented by the radius of k nearest neighbours \(r_{\text{kNN,3D}}\)), and local point density \(\rho_{3D}\). \[H=|Z| \tag{1}\]  
\[\Delta H_{\text{kNN,3D}} \tag{2}\]  
\[\sigma H_{\text{kNN,3D}} \tag{3}\]  
\[r_{\text{kNN,3D}}=max||D_{\text{kNN,3D}}|| \tag{4}\]  
\[\rho_{\text{3D}}=\frac{k+1}{\frac{4}{3}\pi r_{\text{kNN,3D}}^{3}} \tag{5}\]  
where \(D_{\text{kNN,3D}}\) is the distance between the query point and its neighbours. In the other part of the feature extraction, eigenvalue-based features were generated after the neighbourhood was determined. In order to compute these features for a point in the 3D point cloud, the covariance matrix (denoted as **Cov**) needed to be calculated.