2.1.

LiDAR-Image Georegistration

The first step of the fusion process is to perform georegistration between LiDAR and imagery data. We georegistered multispectral LiDAR point clouds to an aerial photo using a phase congruency (PC)-based scene abstraction approach.¹⁹ LiDAR 3D points are converted to two-dimensional (2D)-intensity or height images, whichever type better describes the scene based on the visual interpretation. For example, an area of a wastewater treatment plant varies in height more than in intensity, where most elements are inner and outer tanks, curbs, and structures. These elements share the same concrete cover; and thus, a height-based interpolation would be the optimal decision. On the other hand, a residential area typically includes urban features varying in intensity (i.e., grass, asphalt, and land-markings), which advises an interpolation based on LiDAR data’s intensity values. The approach implements the PC filter, which computes the moment of each pixel’s center point, knowing its PC measure in different orientations. The moment value of a point indicates whether it is an edge or a corner point. A predefined threshold range of moment values can be set to identify candidate tie points on both data sets.

Georegistering data acquired at different times with different sensors may result in two dissimilar sets of candidate tie points, impacting the threshold range. In this case, the PC filter’s outputs (moment images) are abstracted by clustering, then detecting common polylines in LiDAR and imagery data. Moment points within a buffer around these detected edges are considered the candidate tie points. Alternatively, an additional filter can be fused with the PC filter to isolate candidate tie points inputted to the shape context descriptor model to be matched in pairs of final tie points. Finally, a least-squares adjustment estimates the transformation parameters of empirical registration models.

This registration is generic, as it was found to accommodate different urban morphologies, is no longer limited to traditional linear control primitives, and does not require the onboard acquisition of both data simultaneously.¹⁹

2.2.

Color-Based Segmentation

After data registration, LiDAR point clouds are expressed by 3D coordinates as well as their spectral characteristics. The spectral characteristics include the ones originally captured by a LiDAR sensor during the same flight mission and those inherited from aerial images. The color-based segmentation algorithm proposed by Zhan et al.¹⁸ is applied to segment LiDAR point clouds based on their geometric and spectral characteristics. The approach determines the similarity between two points by calculating the geometric and colorimetric distances between one another. It measures colorimetric distances based on the RGB signatures of LiDAR points. It has been applied successfully in previous research work in segmentation-oriented applications.^20,21 The algorithm performs the color-based segmentation in two steps: segment growing and segment merging and refinement.

2.2.1.

Segment growing

In this step, the algorithm assigns each unlabeled LiDAR point to a segment and then marks it as labeled. It constructs three entities: points ($P$) that contains the LiDAR points to be segmented, stack ($ST$), and segments ($S$). In the case of an empty $ST$, which is the default setting when the algorithm starts, the process loops on $P$ until it meets an unlabeled point ($p$). The process appends $p$ to $ST$, creates a new segment ($s$) in $S$, inserts $p$ to $s$, and eventually marks $p$ labeled in $P$. While $ST$ is occupied, it pushes its top point, point of investigation ($pt$) out, and the process searches its neighboring points in $P$ within a 3D distance window. If the neighbors are unlabeled and also radiometrically close to $pt$, the process appends them to $s$ and $ST$. Once $ST$ is clear, densifying $s$ with LiDAR points terminates, and the algorithm looks for another unlabeled point in $P$ to initiate a different segment in $S$. The process continues until $ST$ is empty and all points in $P$ are labeled.

Figure 2 addresses a hypothetical example to illustrate the segment growing step in the color-based segmentation process. Figure 2(a) represents the start run, where all points $P$ are not assigned to segments and yet marked unlabeled, and $ST$ and $S$ are created. $ST$ is by default empty. Hence, the process searches for an unlabeled point ($p_1)$ in $P$, appends it to $ST$, constructs a new segment $s_1$ in $S$, and adds $p_1$ to $s_1$. $ST$’s $pt$ is $p_1$ that becomes the point of investigation [Fig. 2(b)]. The algorithm marks $p_1$ labeled in $P$. $ST$ pushes out its $pt$ ($p_1$), and the process locates its neighbors ($p_2$, $p_3$, and $p_4)$ in $P$ within a predefined 3D distance range [Fig. 2(c)]. The three neighbors are unlabeled, but only $p_2$ meets the radiometric similarity condition. Thus, the algorithm appends $p_2$ to $s_1$ and $ST$ and marks it labeled in $P$ [Fig. 2(d)]. Consequently, the process does not add $p_2$ to a different segment in the future to maintain a one-to-one relation between $P$ and $S$. Since $ST$ is occupied, points inside can contribute to $s_1$’s growth by their neighbor points as long as they are close in color. $ST$ pushes out its $pt$ ($p_2$), which turns to be the point of investigation. The approach determines $p_2$’s neighbors in $P$ ($p_1$, $p_3$, $p_4$, and $p_5$). The latter three points are unlabeled, but none of them meets the colorimetric condition. Therefore, the densification of $s_1$ always ends with an empty $ST$ [Fig. 2(e)].

Fig. 2

Color-based segmentation algorithm applied to LiDAR point clouds: segment growing step. (a) $P$, $ST$, $S$ are created, and $ST$ is empty. (b) Unlabeled $p_1$ in $P$ initiates $s_1$ in $S$, is pushed to $ST$, and becomes $ST$’s $pt$ and point of investigation. (c) $p_1$ is marked labeled in $P$, pushed out by $ST$, and its neighbors are determined in $P$. (d) $p_2$ meets the conditions, is appended to $s_1$, is pushed to $ST$, becomes $ST$’s $pt$, and is marked labeled in $P$. (e) $ST$ pushes out $p_2$ that becomes the point of investigation, and its neighbors determined in $P$. (f) None meets the conditions, $s_1$ is terminated, $ST$ becomes empty, the following unlabeled point in $P$ is $p_3$ that initiates $s_2$ in $S$, $p_3$ is pushed to $ST$, and becomes $ST$’s $pt$ and point of investigation. (g) $p_3$ is marked labeled in $P$, pushed out by $ST$, its neighbors are determined in $P$, qualified neighbors to be examined.

The process loops again on $P$ and locates the following unlabeled point ($p_3$), which initiates a newly created segment $s_2$, joins $ST$, and becomes the point of investigation [Fig. 2(f)]. The algorithm marks $p_3$ as labeled and finds its neighbors in $P$ ($p_1$, $p_2$, $p_4$, and $p_5$). The latter two points are unlabeled; hence, they are the ones for the algorithm to evaluate the colorimetric condition for inclusion into $s_2$ [Fig. 2(g)].

2.2.2.

Segment merging and refinement

The output of the segment growing step is $S$, which contains roughly segmented clusters that need to be merged and refined in this subsequent run. The algorithm builds a merged segments (MS) entity to store lists of homogeneous segments from $S$. The process marks $S$’s segments as labeled if MS already includes those segments. Otherwise, they are marked as unlabeled. The approach first iterates $S$ to find an unlabeled segment ($s$). Then, a new merged segment ($\mathrm{ms}$) is created in MS, with $s$ being a member of $\mathrm{ms}$. Afterward, $s$ is marked as labeled and becomes a segment of the investigation. The process locates its neighboring segments in $S$ within a 3D searching window. $s$’s neighbors are determined by locating the point neighbors of all points within $s$. Every segment to which these point neighbors belong is a neighboring segment. The algorithm calculates the radiometric similarity between $s$ and its unlabeled neighbor segments after averaging the RGB values of each. Suppose that $s$’s neighboring segments are found to be colorimetrically close to $s$, the algorithm appends them to ms and marks them labeled in $S$. The process then continues looping on $S$ and merging segments in MS and terminates when all segments in $S$ are marked labeled and included in MS. Finally, all segments within the same MS are fused. The refinement step looks for MSs with a number of LiDAR points less than a predefined minimum. The process highlights them as MSs of interest, determines their nearest MS in MS within a specific 3D window size, and fuses both in a refined segment (rfs). The refinement continues until the size of all rfs is larger than the predefined minimum. That is when MS eventually turns to be an rfs entity.

Figure 3 graphically explains the segment merging and refinement step. The figure describes a hypothetical scenario where $S$ has five unlabeled segments, and MS is empty before running the step [Fig. 3(a)]. The process loops on $S$, spots an unlabeled segment ($s_1$), marks it labeled in $S$, creates a new MS (${\mathrm{ms}}_1$) in MS, and inserts $s_1$ in ${\mathrm{ms}}_1$ [Fig. 3(b)]. $s_1$ becomes the segment of investigation, for which the process locates its neighbors in $S$ ($s_2$ and $s_3$) [Fig. 3(c)]. Both are unlabeled and meet the colorimetric similarity condition with $s$; hence, the algorithm marks them labeled in $S$ and appends them to ${\mathrm{ms}}_1$ in MS [Fig. 3(d)].

Fig. 3

Color-based segmentation algorithm applied to LiDAR point clouds: segment merging and refinement step. (a) $S$ is from segment growing step, and MS is created. (b) Unlabeled $s_1$ is marked labeled in $S$, ${\mathrm{ms}}_1$ is created in MS, and $s_1$ is inserted in ${\mathrm{ms}}_1$. (c) Segment of investigation is $s_1$, and its neighbors are determined in $S$. (d) $s_2$ and $s_3$ meet the conditions, are marked labeled in $S$, and appended to ${\mathrm{ms}}_1$ in $\mathrm{MS}$. (e) Subsequent unlabeled segment is $s_4$, is marked labeled in $S$, ${\mathrm{ms}}_2$ is created in MS, and $s_4$ is inserted in ${\mathrm{ms}}_2$. (f) Segment of investigation is $s_4$, its neighbors are determined in $S$. (g) $s_5$ meets the conditions, is marked labeled in $S$, and appended to ${\mathrm{ms}}_2$ in MS. (h) Segments are labeled in $S$, merging process terminates, and segments are fused within same $\mathrm{ms}$ in $\mathrm{MS}$. (i) ${\mathrm{ms}}_2$’s size is less than the minimum, it becomes an MS of investigation, its NN (${\mathrm{ms}}_1$) is determined in MS. (j) ${\mathrm{ms}}_2$ is merged to ${\mathrm{ms}}_1$ to form ${\mathrm{rfs}}_1$ in RSF.

The process loops again on $S$, finds the subsequent unlabeled segment ($s_4$), highlights $s_4$ as labeled in $S$, and initiates (${\mathrm{ms}}_2$) in MS with $s_4$ [Fig. 3(e)]. $s_4$ is the segment of investigation, and the process searches its neighbors in $S$ ($s_2$ and $s_5$) [Fig. 3(f)]. Only $s_5$ is unlabeled and meets the radiometric similarity condition with $s_4$. Therefore, the algorithm marks $s_5$ labeled in $S$ and pushes it to ${\mathrm{ms}}_2$ in MS [Fig. 3(g)]. The merging process terminates, since $S$’s elements are entirely labeled and fully included in $\mathrm{MS}$. Segments within same MSs are fused [Fig. 3(h)]. Finally, ${\mathrm{ms}}_2$’s size is less than the minimum, so it becomes an MS of investigation whose nearest neighbor (NN) in $\mathrm{MS}$ is ${\mathrm{ms}}_1$ [Fig. 3(i)]. Eventually, the process merges ${\mathrm{ms}}_2$ to ${\mathrm{ms}}_1$ as a new refined segment (${\mathrm{rfs}}_1$) in $\mathrm{RSF}$ [Fig. 3(j)].

2.3.

Eigenvalue-Based Geometric Features Determination

Sanderson graphically explains the geometric conception of eigenvectors and eigenvalues with a series of visual-aided materials on their website.^22,23 An eigenvector of a linear transformation is a nonzero vector on which the only effect of the linear transformation is scaling by a constant number. The value of the scaling constant is the eigenvalue of the eigenvector. Equation (1) describes the above relation as follows:

To solve for $\lambda$ and $\overrightarrow{v}$, Eq. (1) can be expressed as

Since the right-hand side of Eq. (2) is a zero-vector, $(A-\lambda I)$ is a singular transformation with the property

Equation (3) solves the eigenvalue(s) of the linear transformation $A$. Finally, the corresponding eigenvector of each eigenvalue can be determined using Eq. (2) upon substitution of $\lambda$ with each solved eigenvalue.

Figure 4 shows a numerical example to geometrically describe how eigenvectors and their corresponding eigenvalues are calculated for a transformation matrix $A=\left(\begin{array}{cc}3& 1\\ 0& 2\end{array}\right)$ between 2D spaces. Figure 4(a) assumes the coordinates of the unit vectors $\widehat{i}$ and $\widehat{j}$ in the input space are (1, 0) and (0, 1), respectively. Hence, $(A)$ transforms both to (3, 0) and (1, 2), respectively, which represent the coordinates of the output space’s base vectors with respect to the input space’s grid shown in the background [Fig. 4(b)]. The procedure as mentioned earlier leads to a quadratic polynomial in $\lambda$ that has two solutions: ${\lambda }_1=2$ and ${\lambda }_2=3$, meaning that there are two vectors in the input space that are only scaled by constants upon the transformation. Substituting in Eq. (2) with both eigenvalues determines the corresponding eigenvectors as expressed by $x+y=0$ and $y=0$, respectively [Fig. 4(c)]. A vector that lies on either eigenaxis in the input space does not alter that axis in the output space; the vector is just scaled by its corresponding eigenvalue (2 or 3) as shown in Fig. 4(d).

Fig. 4

Geometric interpretation of eigenvectors and eigenvalues in 2D spaces–numerical example. (a) Coordinate system of input space. (b) Coordinate system of output space. (c) Eigenvectors on input space. (d) Eigenvectors and corresponding eigenvalues on output space.

The geometry of LiDAR points in 3D is derived from the variance of their coordinates: $x$, $y$, and $z$. One way to analyze LiDAR points is to perform a coordinate transformation such that the transformed coordinates $X$, $Y$, $Z$ optimally reveal the variances of the original coordinates. This is achieved when the covariance matrix expressing the variance/covariance between the 3D coordinate records is diagonal, meaning that the diagonal values are the variances of the transformed coordinates $X$, $Y$, and $Z$ that have a zero or minimal correlation between each other. The diagonal covariance matrix in the transformed space ($X$, $Y$, $Z$) can be obtained by eigendecomposition of the original covariance matrix, with each dimension and scaling of the transformed space being an eigenvector and its corresponding eigenvalue of the original covariance matrix. Moreover, the resulted eigenvectors are always orthogonal due to the orthogonality of the transformation (covariance) matrix itself, which is a rotation matrix in this case. The corresponding eigenvalues form indices that better expose the geometry of LiDAR point clouds in a 3D space.^24,25

2.4.

Dimensionality Reduction of Feature Space

It is preferable to diminish the set of input variables (features) of the data being analyzed while developing a predictive model to decrease the computational cost.²⁶ In some cases where the space volume is too large, the data records may not be representative, causing fitting problems that reduce the model performance, in a phenomenon known as the curse of dimensionality.²⁷ Feature selection techniques can be supervised by eliminating statistically irrelevant variables with a weak relationship to the target variable the model attempts to predict. On the other hand, unsupervised feature selection methods drop redundant variables based on statistical measures (i.e., correlation), independent of the target variable.

Even though feature selection and dimensionality reduction techniques attempt a fewer input space to a predictive model, the latter fundamentally differs from the former. Dimensionality reduction methods project the data into a new space of fewer dimensions, resulting in transformed input features²⁶ which are called components in the principal component analysis (PCA) method that we applied in our study to reduce the data dimensionality.

The PCA linearly projects a feature space into a subspace that still preserves the essence of the original data.²⁸ It looks for descriptive features of high variance values revealed by the eigendecomposition of the features covariance matrix. Then, it projects the input feature space into another space constructed by that chosen subset of features. Brownlee²⁸ delineates the process in the following steps:

2.5.

Classification of LiDAR Data

We isolated a portion of the segmented LiDAR data and divided it into training, validation, and testing subsets before employing the multilayer perceptron (MLP) neural network classifier. Table 1 summarizes the main characteristics of the classification we carried out in this work.

Table 1

Description of performed classification.

ANNs are machine learning algorithms that are inspired by observations of how the brain functions with its constituent structures.²⁹ An MLP neural network consists of an input layer and some hidden layers. The last layer is the output layer that provides the final predictions. Each layer is a row of neurons (nodes), where each neuron has weighted inputs, bias, and output to the next layer, representing a perceptron.³⁰ An optimization algorithm weights a neuron’s inputs, and a bias is accumulated to the weighted summation of the inputs to determine their signal strength.³¹ The ultimate goal is setting those weights and biases to particular values so the overall classification error is minimal.³² An activation function is applied to an output signal strength by intensifying or diminishing its value depending on its magnitude. Consequently, outputs of large magnitudes propagate further and contribute to the final predictions more than those of lower magnitudes.³²

Training data have to be numerical in a multiclass classification application. The dimension of the input layer is set to the feature space of the training data, and the output layer has as many nodes as the number of noticed classes. The optimization algorithm initially assigns random weights to each node’s inputs, and their signal strengths are determined after adding a bias to the inputs’ weighted summation.³⁰ The activation function filters signal strengths so only those of high magnitudes propagate to the final predictions. A loss (cost) function measures the model’s classification error, represented by comparing the predictions to their corresponding ground truth data. The model training process iterates over the training dataset for a preset number of times (epochs), or until a condition that signifies a cost function’s minimum is numerically achieved.³¹ The model updates the weights at each epoch, where all training records participate. In addition, internal weights updates occur at each batch or subset of the training data in the case of batch training.³³

3.1.

Study Area and LiDAR Dataset

We used a multispectral LiDAR point cloud captured by the airborne Optech Titan sensor in 2015 to test the proposed approach. The sensor collected LiDAR data using three laser channels: $C_1$, $C_2$, and $C_3$ that represent the 532-, 1064-, and 1550-nm wavelengths of the electromagnetic spectrum, respectively. The point cloud contains 1,976,164 3D points spaced at 0.13 m, produced and projected by the OptechLMS software to the North American Datum (NAD) 1983 Universal Transverse Mercator (UTM) coordinate system (zone 17N). The dataset covers a $\mathrm{33,000}\hspace{0.17em}{\mathrm{m}}^2$ residential area in Rouge within Scarborough, east of Toronto. Megahed et al.³⁴ georegistered the points to a very high-resolution orthophoto acquired in 2014 with a spatial resolution of 0.2 m and covers the same study zone in R, G, B, and NIR bands. They later corrected the overparameterization problem in empirical registration models.³⁵ Figure 5 visualizes the LiDAR dataset before and after the georegistration.

Fig. 5

Visualization of LiDAR data after registration to multispectral very-high-resolution aerial image. (a) Intensity–$C_1$ channel. (b) Intensity–$C_2$ channel. (c) Intensity–$C_3$ channel. (d) Height. (e) R, G, B bands visualized in RGB colors. (f) NIR, B, and G bands visualized in RGB colors.

We used the “lasground-new” tool³⁶ within the collection of LAStools software package to separate ground from nonground points. The tool applies the progressive triangulated irregular network (TIN) densification approach.³⁷ It filtered the LiDAR data into 857,738 and 1,118,426 ground and nonground points, respectively.

We observed four nonground and another four ground classes, as follows:

3.2.

Color-Based Segmentation of LiDAR Data

Figure 6 schematically explains how the segmentation algorithm was applied in the study. It summarizes what Secs. 2.2.1 and 2.2.2 explain in a single chart, with a few alternations that are further illustrated at the end of this section. The classic $k$-NN method is commonly utilized to determine a point’s neighbors. It requires a predetermined $k$ value, representing a chosen number of neighbors for the query point, to initialize the process. First, the algorithm calculates the distances between the point and the remaining data samples. Then, it creates a collection where the indices of those samples are stored with their distances from the query point, sorted in ascending order. Finally, the method selects the first $k$ records from the constructed collection as the point’s $k$-NN.³⁸ Despite the simplicity of the $k$-NN, such a brute-force search is structureless; consequently it is computationally expensive when processing multidimensional data with large sizes.^38,39 Hence, the structure-based $k$-dimensional tree ($k\text{-}d$ tree) method was applied in this study.

Fig. 6

Application of color-based segmentation algorithm with some modifications.

A $k\text{-}d$ tree is a binary tree where each nonterminal node divides the data into two portions depending on a record’s position from a $k\text{-}d$ hyperplane (splitting partition).⁴⁰ Each nonterminal node depicts a different partition. Each level of nonterminal nodes alternates sequentially among the $d$ dimensions in splitting the records. The search starts at the root node and goes down the tree, turning left or right at each nonterminal node based on the query point’s value compared with the threshold value at the split dimension.⁴¹ The search continues until it reaches the terminal node that contains a maximum number of points at which the algorithm switches over to brute-force (leaf size).⁴² However, these points do not represent the final set of NNs for the query point if their extent (centered at the query point) intersects with other hyperplanes. In this case, potential NNs may exist on those sides of the tree where they need to be searched.⁴⁰

A $k\text{-}d$ tree is a data structure based on hierarchical spatial decompositions. Each nonterminal node is associated with a $d$-dimensional cell and the subset of records within this cell. The fundamental design issue is the choice of the splitting hyperplane. The standard split method uses the data distribution by determining hyperplanes orthogonal to the median coordinate along which the points have the most significant spread. However, it generates elongated cells in the case of clustered data, resulting in longer searching times. On the other hand, the midpoint split method uses the cell’s shape by dividing it through its midpoint using a hyperplane orthogonal to the longest side, creating tangibly less elongated cells. Nevertheless, many of the resulted cells can be empty, affecting tree sizes and processing times, especially when dealing with large high-dimensional data. The sliding midpoint method partitions data as the midpoint split method does. However, when it produces empty cells with no data records, it shifts the splitting plane toward data location until the plane touches the first record. This performance overcomes generating sequences of skinny or repetitive blank cells as in the standard and midpoint split methods, respectively.⁴³

We used the “cKDTree” function⁴² embedded in the “scipy.spatial” library in this work. The function constructs a $k\text{-}d$ tree for quick NN lookup. We kept its default values: the sliding midpoint partitioning technique and bucket size of 16 records. Meaning that a node turns to a terminal node (leaf) if 16 points or less are associated with it. Otherwise, the algorithm continues partitioning the data.⁴³

The color-based segmentation algorithm calculates the radiometric distance (RD) between two LiDAR points ($p_1$ and $p_2$) using the Euclidean norm; the square root of the sum of the squares of the differences between the R, G, and B values of each point. We added the NIR figures as follows:

We kept the below threshold values as applied in Ref. 18:

However, we made the following modifications to the algorithm to fit the research objectives:

3.3.

Feature Space Construction

A 10D point features vector represents each point in the LiDAR dataset as below:

These point features fundamentally build the following radiometric and geometric features per LiDAR segment.

3.4.

Feature Space Dimensionality Reduction Using PCA

We used the “PCA” class in the “decomposition” module embedded in the “scikit-learn” library in Python to implement the PCA.⁴⁷ The class implicitly applies the five steps previously mentioned in Sec. 2.4. We created a PCA model and set the number of components to the exact dimensions of the original data. Then, we fitted the model to the input feature space. Afterward, we targeted the components with the highest accumulating variances as the most significant PCs, by which we recreated and refitted the model and finally identified the subspace dimensions. Lastly, we transformed (projected) the original data to the determined output space.

There is no way to name the most significant features within a dimensionality reduction technique; however, the PCA enables recognizing the most meaningful features contributing to each PC’s creation. The “components” attribute in a PCA model is a matrix whose rows resemble the PCs, and columns mirror the features in the input space. Each value in the components matrix represents the weight a feature participated with while producing a PC. By sorting each row in descending order, we could reflect on the most critical features that constructed each PC. We carried out the dimensionality reduction for the ground and nonground subsets separately.

3.5.

LiDAR Data Classification Using MLP Neural Networks

We collected data for classification model training, validation, and testing manually on ERDAS IMAGINE by digitizing polygons of points for each observed class. We acquired 13,953 segments, 80% of which contributed to training and validating the MLP neural network classification model, with a ratio of 80% to 20%, respectively. We used the remaining 20% to test the model (test 1). To ensure a robust model consistency, we added second testing (test 2) using a set of 5337 points collected in the same way. Figure 7 shows the percentage breakdown of the obtained data. It is worth mentioning that the validation and both testing data are three different sets and did not participate in training the model. Also, we maintained a class representation in these four splits proportional to its size in the collected data. Equivalent class representation ensures good model learning and unbiased evaluation. Figure 8 shows the categorization of the eight classes within the collected data.

Fig. 7

Percentage breakdown of collected data: training, validation, and testing.

Fig. 8

Percentage/numerical breakdown of collected data: eight classes. (a) Collected segments (splitted for training, validation, and test 1). (b) Collected points (used in test 2).

Since the class feature in the training, validation, and testing samples belongs to categorical data as a multiclass prediction problem, we converted the class feature to numerical values in two steps. First, we applied label encoding to assign an integer figure to each category (class). Then, we implemented one-hot encoding to give those figures a binary representation to eliminate an ordinal association among the classes.⁴⁸

We used the “Keras” Application Programming Interface (API) in Python⁴⁹ to create an MLP neural network model of four layers. The input layer’s dimension was set to the number of features (37 in total if we entirely consider the radiometric and geometric features). At the same time, the two hidden layers and the output layer had 16, 12, and 4 neurons, respectively. We ran the classification of ground and nonground segments independently; hence, the four neurons in the output layer represent the number of classes observed in both sets. The model was compiled using the Adam optimizer⁵⁰ and the categorical cross-entropy loss function.⁵¹ We applied a traditional feedforward network, where forward processing carried out an upward activation of the neurons until the final output. The rectified linear unit function⁵² activated the input and hidden layers, whereas the softmax function activated the output layer.⁵² Softmax is typically used in multiclass prediction applications, as it provides the probability membership of each segment to belong to each of the output classes. We assigned the class of the highest probability to each segment. The loss function computed the error and backpropagated it, while the optimizer updated the weights according to their contribution to the error. The error backpropagation was iterated over 100 epochs and 50 batches per epoch until the model arrived at a set of weights minimizing the prediction error.³⁰

We examined the following 10 scenarios, each of which trains an MLP neural network model on a bundle of the 37 features, as follows:

Scenarios (1) to (7) thoroughly study the capabilities of the multispectral LiDAR features in urban classification when combined with the height values in seven different ways. Scenario 8 reveals the effect of including R, G, B, and NIR from the aerial image on classifying the scene, whereas scenario 9 tests the hypothesis of combining additional calculated radiometric and geometric features on enhancing the classification results. Finally, scenario 10 attempts a lower dimensionality represented by the most significant PCs instead of the whole input space in the preceding scenario. To ensure a consistent assessment of the 10 scenarios, the training, validation, and testing datasets were the same for the 10 MLP classification models.

3.6.

Classification Assessment

We validated the 10 classification models using different classification metrics and resampling techniques. We constructed the accuracy matrix in test 1 and test 2 for each MLP model to calculate the accuracy, precision, recall, and $F_1$-score metrics. They are explained below for a binary classification that deals with two classes: positive and negative (Fig. 9). We accommodated them accordingly in the calculations to fit a multiclass prediction problem:

Fig. 9

Confusion matrix of binary classification of positive and negative classes.

$\mathrm{TN}$ is the true negative, reflecting the number of negative samples that the model correctly predicts as a negative class, $\mathrm{FN}$ is the false negative, reflecting the number of positive samples that the model incorrectly predicts as a negative class, and $\mathrm{FP}$ is the false positive, reflecting the number of negative samples that the model incorrectly predicts as a positive class,

We trained, validated, and tested the 10 classification models using the same training, validation, and testing sets, respectively, which are illustrated in Sec. 3.5 to maintain an unbiased comparison. In addition, we applied the $k$-fold cross-validation as a different resampling technique. It divides the training records into a $k$ number of folds, where one fold participates as a testing set and the rest contribute to training the model. The algorithm runs $k$ times, and at each turn, a different fold takes part as the held-out testing set. In this way, the entire samples contribute to the fitting and evaluation processes, ensuring robust assessment figures. The $k$-fold cross-validation technique results in a $k$ number of MLP models, whose accuracies are averaged and their standard deviations are calculated.^30,53 In this study, we set $k$ to 10, used the segmented training dataset (Fig. 7), constructed stratified folds to guarantee all classes in training and validation folds are represented proportionally to their size in the training dataset, and repeated the algorithm 100 times.

4.1.

Dimensionality Compression of Feature Space by PCA

The segment growing step in the color-based segmentation process resulted in 38,930 and 92,489 rough segments for the ground and nonground LiDAR points. The merging and refinement step reduced the number of segments to 21,853 and 36,138, respectively, which is 57,991 total segments for the entire point cloud.

After computing the 37D feature space for each segment, we commenced the PCA with 37 PCs, the size of the input space. By sorting the variances in descending order, we calculated the accumulated variance as a preparatory step toward the definite number of components to consider. Figure 10 shows that close to 100% cumulative variance is achieved by recognizing only nine and seven PCs in the analysis of nonground and ground data, respectively. Hence, we reperformed the PCA considering these most significant components and projected the 37D input feature space of the nonground and ground LiDAR points into a 9D and 7D output space, respectively.

Fig. 10

PCs’ cumulative variances. (a) Nonground LiDAR data. (b) Ground LiDAR data.

Figure 11 reveals the contribution of each feature in the 37D input space to creating the most significant PCs after normalizing the features’ weights to range from 0 to 1, indicating no and full contribution, respectively. The most significant PCs of both ground and nonground data primarily consist of combinations of features of direct LiDAR and imagery spectral intensities; ${\mathrm{I}}_1$, ${\mathrm{I}}_2$, ${\mathrm{I}}_3$, R, G, and NIR, with a notable contribution from EVI. As height is an expected feature to distinguish between nonground objects, it is not surprising that it also appears to dominate the most significant PCs of the nonground points. This initial overview highlights anticipation of more influential spectral features than the derived geometric characteristics when the data are classified.

Fig. 11

The most significant PCs’ most contributing features. (a) Nonground LiDAR data. (b) Ground LiDAR data.

4.2.

Classification of LiDAR Data Using MLP Neural Networks

An MLP neural network is a stochastic machine learning algorithm whose decisions vary randomly during the learning process. It uses random initial weights and random shuffle of samples (batches) at each epoch to help the model seek better solutions by avoiding local or deceptive optima. Consequently, each time an MLP neural network algorithm runs, it creates a different model, provides different predictions, and produces a different accuracy. These variations occur even when the algorithm uses the same training data each time it runs.⁵⁴ Hence, it is essential to test the MLP models on different datasets to ensure compatible results.

Figure 12 shows the overall mapping accuracies of the 10 scenarios using the validation, test 1 and test 2 splits (previously addressed in Sec. 3.5), in addition to the $k$-fold approach. In the $k$-fold approach, different nonoverlapping subsets of the training data, which summed up to the entirety of the training set, participated in the validation in rotation. The fourth bar in each scenario represents the mean accuracy of the 10 folds in the 100 repetitions when each fold acted as the held-out testing set. The values above describe the standard deviation of the 1000 accuracy values.

Fig. 12

Assessment of LiDAR data’s multilevel object-based classification: validation, testing, and $k$-fold results.

In each scenario, the four evaluations are close to each other, indicating consistent MLP models that are reliable to apply on unseen data for predictions. However, the first and third models show relatively lower accuracy figures when verified by the test 2 set. This decrease is probably the result of insufficient input feature spaces in both scenarios that allowed misclassified segments, of which some points of test 2 are accidentally part. Test 2 is a set of points, not segments as the rest of the assessment splits are, so the misclassification effect is more pronounced. The evaluation on test 2 shows consistency with the other three assessment splits in the remaining scenarios whose feature spaces are larger with lower standard deviations, which supports this argument. The remaining part of this section discusses the classification results based on the models’ evaluation using the test 1 dataset.

Figure 13 shows the classification results using the test 1 data split. A general glance at the overall mapping accuracy shows its gradual increase by including more features in the input space of the classification. However, there is a significant leap in the nonground accuracy compared with the ground one highlighted in the first three scenarios. This notable increase results from the height being a discriminative feature in the classification of nonground data, whose accuracy is around 90% when combined with a single LiDAR channel. On the contrary, the height range of the ground data is 4 cm, which does not provide room for a better classification when combined with a LiDAR spectral channel, leading to a ground classification accuracy of around 60%.

Fig. 13

Accuracy of LiDAR data’s multilevel object-based classification (test 1 set): ground, nonground, and overall accuracies.

When comparing the first three scenarios, the nonground accuracies of scenario 1 and scenario 2 are almost the same ($\approx 90\%$); however, scenario 3’s is relatively larger (92.37%) because of its higher capability in detecting vehicles. On the other hand, the ground accuracies of scenario 1 and scenario 3 are close ($\approx 61\%$); however, scenario 2’s ground accuracy is relatively larger (64.40%) because of its higher capability in detecting light asphalt. Nevertheless, these variances slightly affect the three scenarios in total; i.e., 79.92%, 80.89%, and 81.68%, respectively.

The following three scenarios reveal the impact of alternating two of the three LiDAR channels. Adding a second LiDAR spectrum booms the overall accuracy to 87.47%, 90.41%, and 91.30% in scenario 4, scenario 5, and scenario 6, respectively, as it relatively enhances the predictions of the vehicles, dark, and light asphalts. The increase of the overall accuracy is vitally contributed by the ground classification, which rushes to around 80% compared with around 60% in the past three scenarios. Combining a second LiDAR beam compensates partially for the idle height feature in the ground classification as per the first three scenarios, with steady progress in the nonground classification figures. Combining the three LiDAR channels in scenario 7 does not add to the highest accuracies of a dual-channel inclusion provided by scenario 6.

Scenario 8 renders another remarkable development in the classification results. Introducing the aerial photo’s radiometric properties (R, G, B, and NIR) increases the accuracy to above 97%. The added spectral features push the nonground and ground accuracies to 97.59% and 97.33%, respectively, striking the overall accuracy of the scenario to 97.49%. This increase results from a continuous improvement in the vehicle, dark, and light asphalt classes.

By accounting for the entire 37D feature space in scenario 9, the nonground classification keeps growing to reach an accuracy of 98.74%. The height-derived geometric feature space allows for better vehicles predictions, raising the accuracy of the nonground classification. However, the full-feature input space slightly affects the dark and light asphalt classes, lightly decreasing the ground accuracy to 97.12%, making the scenario’s overall accuracy 98.17%, and the highest among the entire 10 scenarios.

The most significant PCs in scenario 10 insignificantly lower the nonground classification to 97.91% due to a decrease in the vehicles class accuracy, which is intuitively expected when only a subset of the components participates in the classification. However, the selected components yet contain the most distinguishing characteristics. Unlike the nonground classification, scenario 10 slightly enhances the ground classification to 97.84% due to an increase in the dark and light asphalt classifications. Despite the nontangible improvement, it suggests that the original input space may include a feature or more with a negative impact on the classification model that is eliminated in the PC space, enhancing the results in consequence. Scenario-10 ends with an overall accuracy of 97.89%.

Figure 14 digs more into the classification results by visualizing the per-class accuracies of the different scenarios represented by the $F_1$-score. If a class’s $F_1$-score is zero, it means either its precision or recall is zero. Height and $C_1$ LiDAR channel in scenario 1 can reasonably differentiate buildings, high vegetation, and low vegetation with 95.16%, 80.95%, and 87.61%, respectively, justifying the scenario’s high nonground accuracy (Fig. 13). However, the two features show buildings/high-vegetation and vehicles/low-vegetation misclassification problems.

Fig. 14

Per-class accuracies of LiDAR data’s multilevel object-based classification: $F_1$-score (test 1 set).

Because of the different data acquisition times (LiDAR and imagery data were obtained independently in 2015 and 2014, respectively), the model misclassifies many high vegetation segments as buildings. Some high-vegetation LiDAR points resemble asphalt locations on the aerial image, where trees were not planted. The model also misclassifies some segments of buildings as high vegetation due to orthorectification problems. These problems are attributed to the building facades appearing in the aerial image; consequently, LiDAR segments on roof edges inherited the spectra of building interfaces and surrounding surfaces usually planted in residential areas.

More severely, scenario-1 shows poor vehicle classification results (47.76%) as it mixes a substantial part of the class with low vegetation. It is unlikely for a vehicle to keep the exact location in a study scene, particularly when captured by two different sensors on different dates. Therefore, they inherited the corresponding spectra of the ground surfaces where they usually park (i.e., grassy sidewalks). Hence, the abovementioned pairs of misclassified classes are radiometrically and geometrically indistinguishable, with the height information being the solely geometric feature in scenario 1.

On the other hand, scenario 1 recognizes sidewalks and grass with an accuracy of 88.29% and 81.85%. Nevertheless, it reports dark-asphalt/light-asphalt/grass misclassification that plunges dark asphalt to 31.14% and completely misses light asphalt. Scenario 2 introduces $C_2$ LiDAR channel instead of $C_1$. It increases the high-vegetation accuracy (85.08%), almost does not change the accuracy of buildings (95.89%), decreases the accuracy of low vegetation (85.60%), and completely drops the accuracy of vehicles due to the previously mentioned misclassifications, which eventually does not alter the nonground classification as noticed in Fig. 13. Nonetheless, scenario 2 identifies light asphalt with an accuracy of 66.56% after a drop in scenario 1. This jump enhances the ground accuracy (Fig. 13) despite the misclassification of the entire four ground classes, which misses the sidewalks and dark asphalt, and lowers the grass accuracy to 78.10%.

Scenario 3 tests $C_3$ instead of the other two LiDAR channels: $C_1$ and $C_2$. The model boosts the vehicle classification to 64.15% after a drop in scenario 2, consequently raising the low vegetation to 90.59%. By providing a slight increase to the classification of buildings (96.40%), scenario 3 increases the nonground classification accuracy compared with scenario 1’s and scenario 2’s (Fig. 13). It also records a hit in classifying sidewalks (98.10%), compensating for missing both asphalts as grass, lowering the grass accuracy to 72.78%, in consequence. This hit keeps scenario 3’s ground accuracy close to scenario 1’s.

Alternating the three LiDAR channels in pairs as per scenario 4 to scenario 6 presents a nearly steady increase of all classes, converging in a high accuracy range, from 87% to 98%, except for vehicles and asphalts. The inclusion of dual channels rushes the accuracy of light asphalt to 73.89%, 76.16%, and 76.47% in the three scenarios, respectively. The dark asphalt class accuracy also jumps from 2.84% in scenario 4 to 30.77% and 29.47% in scenario 5 and scenario 6, respectively. On the other hand, the vehicles’ accuracy decreases to 45.33% in scenario 4, then raises to 59.06% and 68.16% in scenario 5 and scenario 6, respectively. These figures give scenario 6 higher ground, nonground, and overall accuracy than the preceding scenarios and are very close to scenario 7’s, donated by including the three LiDAR channels (Fig. 13). Compared with scenario 6, scenario 7 decreases vehicles from 68.16% to 58.99% and increases dark asphalt from 29.47% to 33.33%.

The last three scenarios continue the nearly steady improvement of all classes, yet converging in a higher accuracy range, from 96% to full prediction, except for vehicles and both asphalts. Accumulating the aerial photo’s spectra in scenario 8 tangibly improves the troubling classes to 88.30%, 90.85%, and 95.13% for vehicles, dark, and light asphalt, respectively. The R, G, and B bands lend a hand in discriminating dark and light asphalts even with the naked eye visualization [Fig. 5(e)]. At the same time, the NIR band helps identify the greeny features with 96.98%, 96.62%, and 99.60% accuracy for the high, low vegetation, and grass classes, respectively, which raises the accuracy of the vehicles and buildings (98.80%) in accordance since they are misclassified with low vegetation and high vegetation, respectively.

The ${\mathrm{SF}}_{\mathrm{geom}}$ introduced in scenario 9 better eliminates this misclassification, as the mixed classes share close values of the ${\mathrm{SF}}_{\mathrm{rad}}$ added in the same scenario since they already have similar green characteristics. Scenario-9 notably increases the vehicle accuracy to 92.55% and develops the accuracy of buildings, high vegetation, and low vegetation to 99.55%, 98.58%, and 97.85%, respectively. On the other hand, the scenario suggests one or more confusing features in the ${\mathrm{SF}}_{\mathrm{rad}}$, such as lowering the sidewalks to 97.65% after a full prediction in scenario 8, also decreasing the light asphalt to 94.29%. The PCA eliminates confusing features by definition; therefore, scenario 10 increases the accuracy of dark, light asphalts, sidewalks, and grass to 92.94%, 95.54%, 100.00%, and 99.80%. These figures are not only higher than scenario 9 but also even exceed what scenario 8 achieves for ground classes.

4.3.

Production of Final Urban Map

The bar charts in Fig. 15 summarize the best and the worst classification accuracy results we achieved in this study from scenario and class perspectives. Figure 15(a) displays the most useful versus the most unfavorable scenarios for the land-uses, ground, nonground, and overall classification, based on the results discussed in Sec. 4.2. Scenario 9 produces the highest nonground accuracy (98.74%) as expected after the scenario being the highest-scoring in the nonground classes, except for high vegetation. The class yields 98.94% accuracy in scenario 10, slightly above the 98.58% obtained by scenario 9. Likewise, scenario 10 provides the best ground (97.84%) and ground per-class accuracies, except for grass that reaches 99.90% in scenario 9, insignificantly over the 99.80%, it hits with scenario 10. Therefore, we considered scenario 9's and scenario 10's nonground and ground classifications to produce the final urban map.

Fig. 15

Summary of LiDAR multilevel classification’s best and worst achieved results: $F_1$-score (test 1 set). (a) Scenario perspective. (b) Class perspective.

Figure 15(b) shows the same results but from a different aspect, as it sums up the most versus the least detectable classes for each scenario. Building roofs are the best hits of the majority of the scenarios. As long as a feature space includes the height records, adding a single radiometric feature guarantees $>95\%$ detection (i.e., 95.16% in scenario 1), which can be improved by adding more features. Scenario 3 is the best to efficiently target sidewalks or similar elements (i.e., landmarking at pedestrian crossing intersections) with a $>98\%$ accuracy.

Scenario 2 and scenario 4 are still easy-to-pick options (narrow feature spaces, and thus, fast processing) if a building-accuracy higher than what scenario 1 offers is required. Scenario 4 is also beneficial in grass-oriented applications with a $>92\%$ accuracy (Fig. 14). Scenarios 5 to 7 are good choices if a $>97\%$ building detection is needed or a $>90\%$ green accuracy is attempted ($>90\%$, $>92\%$, $>95\%$ for low, high vegetation, or grass, respectively) [Fig. 14)].

Scenario-8 is a perfect compromise of the entire eight classes. Besides a full grass detection, the inherited R, G, B, and NIR substantially solve the vehicle/low-vegetation and dark/light-asphalt misclassifications emerging in the preceding scenarios. Vehicles are the scenario’s least accurate features, yet, detected with an accuracy of 88.30%. Nevertheless, scenario 9 is optimum for the maximum misclassification elimination between geometrically distinguished classes (i.e., vehicle/low-vegetation and building/high-vegetation). Consequently, the scenario fits urban mapping applications with fine accuracy requirements ($>98\%$). Vehicles are still the scenario’s least accurate features, but 91.39% accurate. Scenario 10 suits large input feature spaces when one lacks a predetermined knowledge about their significance.

We want to emphasize that the results’ discussions, along with the recommended case uses of each scenario, are guidelines for researchers, assuming similar urban objects and encountered data challenges (i.e., orthorectification, shadow, and different acquisition time). Researchers should consider their data structure, observed classes, and application requirements to decide optimal classification scenarios. Our analysis may shed light on even more feature combinations and testing scenarios.

Picking the best scenario depends on the application’s nature and objectives. This study aims to provide the best possible accurate urban mapping that accommodates all accuracies: per-class and overall accuracies. Consequently, we chose scenario 9 and scenario 10 for nonground and ground classifications. The combined scenarios increase the overall final map accuracy to 98.43%, slightly higher than the maximum overall accuracy achieved by scenario 9 alone [98.17%; Fig. 15(a)].

Figure 16 shows the per-class accuracy of the combined-scenario classification. It separates the nonground and ground classes accuracies of scenario 9 and scenario 10, respectively, from Fig. 14 for better visualization of the final map’s per-class figures visualized in a bar chart.

Fig. 16

Produced urban thematic map’s per-class accuracies: $F_1$-score (test 1 set).

Figure 17 reflects on the learning curves of the MLP neural networks of the ground and nonground classifications, which we used to produce the final urban map. Both models show good performance on training and validation samples since both reach minimal losses (errors). However, good learning is revealed here by the loss convergence of the training and validation datasets around a relative value each time the epochs increase. On the contrary, overfitting occurs when the validation curve deviates with higher loss values from the training curve after a convergence. In comparison, underfitting happens when training data always show lower loss values than the validation samples. In this case, the validation curve either declines or levels off with the increase in epochs.⁵⁵

Fig. 17

Produced urban thematic map’s MLP learning curves. (a) Nonground model (scenario 9). (b) Ground model (scenario 10).

Figure 18 shows the classified LiDAR point cloud as the final produced urban thematic map, also used for the qualitative assessment of the classification. The map [Fig. 18(a)] shows the eight classes accurately placed as the quantitative evaluation results suggest, in comparison with the corresponding aerial image [Fig. 18(b)] used in the data registration process. We highlight on the map example locations of five misclassification cases explained as follows:

Fig. 18

Produced urban thematic map. (a) Classified point cloud with examples of misclassification locations: (1) vehicles/low-vegetation; (2) shadow and different acquisition time of data; (3) buildings/high-vegetation; (4) objects with no class representation; and (5) buildings/high-vegetation. (b) Corresponding aerial image: R, G, and B bands visualized in RGB colors.⁵⁶

Nonetheless, our results based on the qualitative and quantitative assessments outperform the nonground and overall accuracies achieved by Megahed et al.,³⁵ who carried out a point-based classification with scenario 8 on the same LiDAR point cloud using the same classifier. This comparison underlines the efficiency of introducing height-derived geometric features in classifying urban objects that vary in height values.