A 3 D Shape Descriptor Based on Contour Clusters for Damaged Roof Detection Using Airborne LiDAR Point Clouds

The rapid and accurate assessment of building damage states using only post-event remote sensing data is critical when performing loss estimation in earthquake emergency response. Damaged roof detection is one of the most efficient methods of assessing building damage. In particular, airborne LiDAR is often used to detect roofs damaged by earthquakes, especially for certain damage types, due to its ability to rapidly acquire accurate 3D information on individual roofs. Earthquake-induced roof damages are categorized into surface damages and structural damages based on the geometry features of the debris and the roof structure. However, recent studies have mainly focused on surface damage; little research has been conducted on structural damage. This paper presents an original 3D shape descriptor of individual roofs for detecting roofs with surface damage and roofs exhibiting structural damage by identifying spatial patterns of compact and regular contours for intact roofs, as well as jagged and irregular contours for damaged roofs. The 3D shape descriptor is extracted from building contours derived from airborne LiDAR point clouds. First, contour clusters are extracted from contours that are generated from a dense DSM of individual buildings derived from point clouds. Second, the shape chaos indexes of contour clusters are computed as the information entropy through a contour shape similarity measurement between two contours in a contour cluster. Finally, the 3D shape descriptor is calculated as the weighted sum of the shape chaos index of each contour cluster corresponding to an individual roof. Damaged roofs are detected solely using the 3D shape descriptor with the maximum entropy threshold. Experiments using post-event airborne LiDAR point clouds of the 2010 Haiti earthquake suggest that the proposed damaged roof detection technique using the proposed 3D shape descriptor can detect both roofs exhibiting surface damage and roofs exhibiting structural damage with a high accuracy.


Introduction
Building collapse is one of the primary causes of heavy human casualties in destructive earthquakes [1].Rapid and reliable damage assessment on the individual building level following earthquakes has become imperative for the optimal utilization of available resources for rescue [2,3].Damage to roofs is an important feature for distinguishing extreme damage states, i.e., collapsed buildings, from lesser damaged or undamaged buildings [4].Therefore, vertical remote sensing, including optical, SAR and LiDAR, represents an efficient tool for rapid damage assessment due to its low cost, high availability, minimal corresponding fieldworks, large coverage, digital processing and quantitative results.The effectiveness of remote sensing has also been proven following earthquakes worldwide [5][6][7][8].Numerous methods have been reported for building damage detection using 2D features, such as gray scale, spectra, texture, edge and morphological features, and amplitude and phase information derived from optical or SAR imagery [8][9][10][11][12][13][14][15][16][17].However, certain damage types (e.g., pancake collapse) cannot be identified using vertical remote sensing data due to the absence of precise height data [18].Airborne laser scanning systems are particularly suitable for damaged roof detection because precise 3D point clouds can be rapidly obtained at all times and under most weather conditions without entering the quake-stricken area [19,20], and the elevation accuracy is higher compared to point clouds derived from vertical optical or SAR imagery [18,21].Change detection using both pre-and post-event remote sensing data is a popular method of acquiring building damage information because detailed pre-event data are invaluable in reconnaissance [18,[22][23][24].However, the major limitation concerning this method is the lack of homogeneous pre-event reference data in many situations [18].The method of damage interpretation using post-event remote sensing data can be applied even in the absence of homogeneous reference data, which is an alternative to rapid damage assessment during earthquake emergency response when pre-event data are limited.
According to the essential clue that damaged buildings, unlike organized manmade patterns on intact buildings, usually manifest themselves as disturbed spatial or spectral patterns [25], various methods have been developed to infer damage patterns from post-event data.The 2D features, including edge, texture and spectra, have been assessed by numerous studies [17,[26][27][28] as important cues for damage detection because damaged regions tend to exhibit disturbed spatial or textural patterns, in contrast to intact buildings [29].On the other hand, 3D features have been found to be useful for identifying specific damage types based on geometric reasoning, as highlighted by some studies [30][31][32][33].Therefore, this paper focuses on severely damaged buildings and proposes a damaged roof detection approach based on the roof's 3D features extracted from only post-event airborne LiDAR point clouds.In the following subsection, we do not provide an exhaustive review of all these methods; instead, we highlight only the 3D-feature-based approaches using only post-event airborne point clouds that are directly relevant to our work in the next subsection.

Damage Types
Prior knowledge about damaged buildings is necessary when performing building damage detection using airborne LiDAR point clouds [34].A building damage catalog, as shown in Figure 1, including typical damage types of buildings, is used to identify buildings exhibiting different damage types based on geometry features such as reductions in volume and height, changes in the inclination of building surfaces and surface structures as well as the size of debris [35].However, extracting geometry features such as reductions in volume and height and changes in the inclination of building surfaces using only post-event airborne LiDAR point clouds is difficult.Therefore, a majority of the building damage detection approaches using only post-event airborne LiDAR point clouds utilize geometry features including the surface structure and debris.This paper categorizes the building damage catalog into two categories based on geometry features of roofs including debris and the roof structure.The first category, named surface damages, contains damaged building with debris surfaces such as multilayer collapse (2), top story pancake collapse (4c, 5, 5c), heap of debris (6, 7a, 7c) and heap of debris with planes (3, 7b).The second category, named structural damages, contains damaged buildings with relatively intact surfaces and includes inclined plane (1), middle or lower story pancake collapse (4a, 4b, 5a, 5b) and inclination (9a).Pancake collapse is a damage type of concern.The top story pancake collapse types, including damage types 4c, 5 and 5c, can be detected using post-event airborne LiDAR point clouds because the roofs are collapsed or damaged, whereas middle or lower story pancake collapses, including damage types 4a, 4b, 5a and 5b, are difficult to detect because the roofs are nearly intact.However, few perfect middle or lower story pancake collapses wherein the building entirely maintains its surface and structure during destructive earthquakes occur; the majority of buildings exhibiting this type of pancake collapse usually exhibit an inclination as well [36].Therefore, this form of pancake collapse can be detected through the structure analysis of inclination.

Building Damage Detection Approaches
The majority of approaches to damage detection using only post-event airborne LiDAR data can be categorized into surface damage detection approaches and structural damage detection approaches, according to the damage categories defined in the above section.
The majority of methods for surface damage detection detect the surface damages of collapsed buildings based on the planarity of the roof surface because airborne LiDAR point clouds are particularly suited to extracting planar roof surfaces [37][38][39].Rehor et al. [40] produced a 2.5D planar Delaunay-based triangulated irregular network (TIN) using planar segments and non-segmented points segmented from roof points and extracted debris triangles from the TIN as the damaged building parts.Rehor et al. [41] compared the random sample consensus (RANSAC) and region-growing algorithms applied to Digital Surface Models (DSMs) for building damage detection and suggested that the region-growing algorithm is more suitable for building damage detection than the RANSAC algorithm.Labiak et al. [42] presented a line-based slope threshold method for evaluating and identifying the damaged points of each roof based on the idea that points in intact roof planes have constant slopes, whereas points in damaged roof surfaces have varying slopes.Segment-based classification methods have been presented to detect collapsed buildings using only post-event airborne LiDAR point clouds based on the assumption that damaged buildings are represented as many small planar segments or unsegmented points, whereas intact buildings are represented as large planar segments [43][44][45].
Structural damage detection approaches detect damaged buildings via structure analysis based on prior knowledge about intact buildings.Shen et al. [31] extracted the geometric axis line of a flat or symmetric roof and identified the inclined roofs based on the assumption that the inclined roof's angle between its geometric axis line and plumb line is greater than an empirical threshold angle.Gerke and Kerle [32] presented a graph-based approach for structural seismic damage assessment based on oblique airborne images.The structural integrity of a building is inferred based on the spatial relation between observable features such as vegetation, façade, intact roofs and destroyed roofs.The relations are represented through a directed graph and are trained by a graph-based learning algorithm.However, this approach is difficult to apply directly to airborne LiDAR point clouds because it is difficult to extract building façades from airborne LiDAR point clouds.Vetrivel et al. [29] developed a gap-based classification method for building structural damage assessment using post-event image-based 3D points.In this method, 3D points of building elements are voxelized based on a pre-defined voxel size; then, gaps are identified as the voxels that are visible from a sufficient number of cameras but that are not occupied by 3D points.Finally, using radiometric features, the gaps caused by damage are detected based on surrounding damage patterns such as spalling or debris.However, it is difficult to extract such evidence from airborne LiDAR point clouds because of their limited radiometric information.As a result, damage-related gaps cannot be reliably classified using airborne LiDAR point clouds.Fernandez Galarreta et al. [4] presented an UAV-based method for urban structural damage assessment using object-based image analysis and semantic reasoning.In this method, the detailed 3D point clouds were generated from multi-view imagery obtained by unmanned aerial vehicles; then, the z component of the normal of a local tangent plane of each point was computed from the co-variance matrix of the neighborhood points and was used to visually assess the D4-D5 damage elements in terms of the European Macroseismic Scale 1998 (EMS-98).Finally, the D1-D3 damage features were extracted via a more detailed façade and roof analysis using object-based image analysis and semantic reasoning.However, the D4-D5 damages were identified not by automatic analysis but by visual analysis.
The above-mentioned methods using airborne LiDAR point clouds all mainly focus on surface damage and pay minimal attention to structural damage.Thus, the spatial relations between the components of a building must be analyzed to infer structural damage [32].The topological relationships between adjacent planar segments can be described and reconstructed using a roof topology in building modeling [46][47][48].However, it is challenging to mathematically describe and explain the topological relationships of complex, damaged roofs based on a scattered point cloud at a low level, including at the point level and segment level, because the topology of random and irregularly damaged roofs is uncertain.
The 3D shape of a complex, damaged roof can be quickly represented by contours, which avoids the problem of topological analysis and reconstruction based on planar segments [49].By analyzing the characteristics of a building's contours, some of their features, such as closed and regular shapes, simple topology and density, are leveraged to extract and reconstruct buildings [49][50][51][52].Our previous studies have shown that it is possible to detect roof damage using a contour-based method [53].Damaged roofs with confusing contours were detected by a shape similarity analysis algorithm applied to building contours derived from airborne LiDAR point clouds.However, the damage feature based on contours was not explicitly described, and the automation of the algorithm was poor.
Focusing on surface and structural damages on roofs, this paper defines a 3D shape descriptor based on shape analysis of the contour clusters of buildings for detecting severely damaged buildings from post-event airborne LiDAR point clouds.The contribution of the paper lies in the presentation of a 3D shape descriptor that provides a comprehensive description of both surface and structure features of roofs based on the shapes and spatial relations of building contours.Compared to other 3D features of roofs, the 3D shape descriptor can more reliably and completely detect damaged roofs.
The remainder of this paper is organized as follows: Section 2 details the key procedures of damaged roof detection, including the data preprocessing, the definition and calculation algorithm of the proposed 3D shape descriptor, and damaged roof detection based on the 3D shape descriptor using a maximum entropy threshold.Section 3 introduces the study area and the data source.Section 4 presents experimental results and discussion about damaged roof detection on the airborne LiDAR dataset of the Haiti earthquake in 2010.The selection of values for algorithmic parameters is also discussed in Section 4. The conclusions are presented in Section 5.

Methodology
This study proposes a method of damaged roof detection for both roofs suffering surface damage and roofs suffering structural damage based on a 3D shape descriptor using post-earthquake airborne LiDAR point clouds.The 3D shape descriptor is used to quantitatively describe both surface and structure features of roofs based on the shapes and spatial relations of a building's contours.The descriptor is a more comprehensive description of 3D shapes of damaged roofs compared to other 3D features.Therefore, damaged roof detection based on the 3D shape descriptor can be more reliably and effectively performed.
There are three key procedures that constitute this method: data preprocessing, feature extraction and damaged roof detection.In data preprocessing, the DSM of each individual building is extracted from post-earthquake airborne LiDAR point clouds with guidance from 2D GIS vector data of building footprints and a digital elevation model (DEM).In feature extraction, the 3D shape descriptor of each building is calculated based on contour clusters generated from the DSM of each individual building.In damaged roof detection, damaged roofs are detected based on the 3D shape descriptor.Figure 2 illustrates the detailed procedures.

Data Preprocessing
With the guidance of the 2D GIS vector data on building footprints, a dense DSM of an individual building is constructed using the airborne LiDAR point cloud and the DEM.The individual building points are extracted from the point cloud according to the building footprint, as observed in Figure 3b,f.A thin plate spline (TPS) interpolation method is used to interpolate the dense DSM with a grid cell size of λ to remove the impact of noise and to create a continuous and smooth surface [54,55], as shown in Figure 3c,g.The ground points are set as having the same height, which is acquired from the DEM.

Feature Extraction
The shapes and spatial relations of building contours can represent the roof's 3D shape feature integrated with both surface and structure features more effectively than can planar segments.A 3D shape descriptor is defined to quantitatively describe the 3D shape feature using the shapes and spatial relations of building contours.Therefore, the 3D shape descriptor can potentially be used to identify surface or structural damages on roofs.

Feature Definition
The 3D shape of a roof's surface and structure can be represented by a contour cluster, which is a set of contours extracted from a contour tree based on a containment relationship.A contour tree is a data structure addressing the description of surface topologies and is built upon the containment relationship between contours [56].A contour cluster is defined as a subtree in the contour tree according to the following three conditions: (1) the subtree's root node has no parent node or has a brother node in the contour tree; (2) the subtree's terminal node has no child node or has at least two child nodes in the contour tree; and (3) the subtree's other node has only one child node in the contour tree.The building's surface regions belonging to the same structure are segmented from the complex building surface by the contour cluster based on a homogeneous spatial topology relationship, and each region is represented by a set of contours in the contour cluster, as shown in Figure 4.The 3D shape difference between damaged and intact roof surface regions can be reliably distinguished using the full 3D geometric relationships among the contour clusters.One of the essential differences between damaged and intact buildings is that the former usually manifest as disturbed spatial patterns such as irregular surfaces and structures, unlike the organized, manmade patterns of the latter [25], as demonstrated in Figure 5.For an intact roof, the organized manmade pattern of the regular surface and structure is represented by a contour cluster containing similar shape contours.For a damaged roof, the disturbed spatial pattern of the irregular surface and structure is represented by a contour cluster containing chaotic shape contours.Therefore, the shape chaos of a contour cluster can be used to describe the irregular 3D shape of a damaged roof's surface and structure and thus represents a potential feature for reliably distinguishing between intact roofs and damaged roofs due to surface or structural damages.According to the concept of information entropy developed by Shannon [57], the shape chaos index is defined to quantitatively describe the shape chaos based on the contour shape similarity measurement.The principle of information entropy is to use uncertainty as a measure to describe the information contained in a source [58].For example, entropy is used to quantitatively measure the spatial information of a map such as metric information, topological information and thematic information [59].In this paper, entropy is used to quantitatively measure the shape chaos of the entire contour cluster based on the shape similarities between two contours.
This can be achieved using the ratio between the number of contours with similar shapes and the total number of contours in the contour cluster as the probability P i used in the definition of the entropy E, as shown in Equation ( 1).Let N be the total amount of contours, clustered into n groups with similar shapes, and let the number of contours in each group be N i .Thus, the probability P i can then be defined as shown in Equation (2).Specifically, if each contour is sufficiently similar, the shape chaos index will equal zero.If each contour is very different, which leads to each group having only one contour, the contour cluster will be clustered into N groups, and all probabilities will be 1/N; therefore, the shape chaos index will reach a maximal value: ln(N).Thus, the shape chaos index C is the entropy E normalized by the maximal value, as shown in Equation ( 3).
E " ´ÿ n 1 P i ln pP i q (1) C " ´řn The 3D shape descriptor is the weighted sum of the shape chaos indexes of contour clusters corresponding to a single building, which describes the chaotic 3D shape of the entire building surface.The 3D shape descriptor is defined as shown in Equation (4).

S "
where m is the number of contour clusters corresponding to a single building, Q i is the weight of each contour cluster, and C i is the shape chaos index of each contour cluster.In this paper, Q i can be determined using the ratio between the area of the building surface region segmented by the contour cluster and the total area of the building.Let A be the total area, and let A i be the area of each region.Such a weight can then be defined as shown in Equation ( 5).Because the shape chaos index and the weight are both normalized, the 3D shape descriptor is also normalized.
The damaged roof can be explicitly described by the 3D shape descriptor based on contour clusters.The 3D shape descriptor measures the chaotic degree of the entire roof surface's 3D shape based on the shape chaos index.According to the Shannon entropy [57], the greater the difference between the contour shapes within the contour cluster, the larger the shape chaos index of the contour cluster, which means that the 3D shape of the surface region segmented using the contour cluster is more irregular.If the area of the irregular surface region is greater, the 3D shape descriptor will be larger.This means that the probability that the roof was damaged is higher.Accordingly, we assume that damaged roofs can be distinguished from intact roofs by a threshold δ of the 3D shape descriptor.Furthermore, the damaged roof is defined using the 3D shape descriptor, as shown in Equation (6).

Feature Extraction Algorithm
This paper proposes a feature extraction algorithm for a 3D shape descriptor based on the shape analysis of contour clusters.There are three key procedures in the feature extraction algorithm: the contour cluster extraction, the shape chaos index calculation and the 3D shape descriptor calculation.Contour clusters are extracted from the contour tree, which is built using relationships based on the containment method for contours generated from the DSM.The shape chaos index of the contour cluster is calculated based on the information entropy of shape similarities, which are computed based on the shape similarity measurements among contours.The 3D shape descriptor is calculated based on the shape chaos indexes.Figure 6 shows the complete workflow of the feature extraction algorithm.(1) Contour Cluster Extraction Contour clusters are extracted from the contour tree based on the three conditions discussed in Section 2.2.1.Dense contours with a contour interval of ε are generated from a dense DSM of an individual-building-based grid-tracking method that interpolates points at given heights between grid cells and connects them sequentially [60,61], as shown Figure 3d,h.The contour tree is built using relationships based on the containment method [62,63].The containment method suffers from the limitation that all contours are considered closed, and the unclosed contours are thus initially excluded.In a contour tree, each node represents a different contour, and every node may have a list of descendants.If one contour is contained by another, then that contour is a descendant.The containment method begins from the root node with the lowest elevation and recursively creates the contour tree in a depth-first manner.
(2) Shape Chaos Index Calculation The shape chaos index of the contour cluster is used to quantitate the chaotic spatial pattern as the information entropy of shape similarities among contours.The shape similarity measurement is a crucial step in calculating the shape chaos index.The shape similarity is measured using the Euclidean distance between the normalized Fourier descriptors (FDs) of two contours.FDs are typically used in image retrieval and pattern recognition because of their insensitivity to geometric translation, rotation, and scaling [64][65][66][67].We employ the normalized Fourier descriptor (nFD) to measure the approximate shape similarity between two contours based on the assumption that the shape differences among the intact contour clusters are only produced by geometric translations and scaling, whereas the damaged contour clusters do not follow this assumption.Thus, the shape similarities based on nFDs among the intact contour clusters will remain small, and the shape similarities among the damaged contour clusters will be random.This results in the shape chaos indexes of contour clusters of intact roof surface regions having small values, whereas such values will be large in damaged roof surface regions.
Suppose a contour is composed of M points as a discrete complex function; then, the function can be transformed into the frequency domain without any loss of information by the Discrete Fourier Transformation (DFT), defined as a(k) in Equation (7).The coefficients a(k) are FDs that represent the discrete contour of a contour line in the Fourier domain [64].
The Fourier coefficients a(k) are normalized to make them invariant to the translation, rotation, and scaling of contours [67].The normalized FDs are defined as in Equation (8).Translation invariance can be achieved by omitting the Fourier coefficient a(0) and using the other Fourier coefficients because translation affects only the first Fourier coefficient a(0), and the other FDs retain their values [64].The FDs are made invariant against rotation by taking the magnitude of each Fourier coefficient, and they are made scale invariant by dividing all Fourier coefficients by the magnitude of a(1) [64,67].The normalized FD na(k), defined in Equation ( 8), also omits the second Fourier coefficient because it is a constant value of 1.
The Euclidean distance is based on the normalized FDs of two contours following Equation ( 9).
s " where na α and na β are the normalized FDs of contour α and contour β, respectively, and L is the number of normalized Fourier coefficients.In most cases, these two contours contain a different number of points, and their normalized Fourier coefficients are not same; therefore, the Euclidean distance cannot be directly calculated.To solve this problem, both contours are first sampled such that they have the same number of points before FDs are applied to the two contours [67].
The shape chaos index is calculated based on an entropy definition whereby the probability is defined as the ratio between the number of contours with similar shapes and the total number of contours in the contour cluster, as detailed in Section 2.2.1.The first step is to cluster contours into several groups following the clustering rule that the shape similarity between two contours in the group must be less than a threshold ω.Based on the results of contour clustering, the probabilities of the shape distributions are computed, and the shape chaos index is calculated following Equation (3).
(3) 3D Shape Descriptor Calculation The 3D shape descriptor calculation is dependent on the shape chaos indexes of contour clusters and their area weights.The area weight is defined as the ratio between the area of the roof surface region segmented by the contour cluster and the total area of the roof.Thus, it is crucial to estimate the area of the roof surface region for the weight calculation.The contour cluster is used to estimate the area of the corresponding roof surface region.Based on the containment relationship between contour clusters, there are two types of contour clusters.The first type of contour cluster does not contain any other contour cluster and has the same area as its outermost contour.The other type of contour cluster fully contains some contour clusters.The area of the cluster is the area of the outermost contour reduced by the areas of the outermost contours in contour clusters contained by the contour cluster.The total area of the roof's surface is the sum of the areas of all contour clusters.Finally, the weight is computed following Equation (5).

Damaged Roof Detection
The damaged roof detection problem is transformed into a roof classification problem by using a maximum entropy threshold for the 3D shape descriptor.The maximum entropy threshold is implemented to select the optimal threshold δ from a group of buildings that includes intact and damaged buildings.As an optimal criterion, the maximum entropy was initially used for image thresholding by Pun [68] and was later corrected and improved by Kapur [69].The maximum entropy threshold considers an image histogram as a probability distribution, and later, one obtains an optimal threshold value by maximizing the upper bound of the a posteriori entropy [58,70].This paper adopts Kapur's method based on the Shannon entropy to select an optimal threshold δ that yields the maximum entropy between damaged and intact buildings.
Suppose that all roofs are classified as intact or damaged roofs using a threshold of the 3D shape descriptor, which is denoted as t.Thus, the numbers of intact and damaged roofs are determined by the threshold t and are defined as functions of t, as shown in Equations ( 10) and (11).
where N I is the number of intact roofs and N D is the number of damaged roofs.Intact and damaged roof histograms are built based on the 3D shape descriptor with a bin width ϕ.Because the 3D shape descriptor is normalized, the number of bins is denoted as B, which is determined using Equation ( 12).The frequency of each bin of the intact roof histogram is FI i , as given in Equation ( 13), and the frequency of each bin of the damaged roof histogram is FD i , as given in Equation ( 14).
B " The entropy of the intact roofs is E I , as given in Equation ( 15), and the entropy of the damaged roofs is E D , as given in Equation ( 16).The a posteriori entropy is defined as a function of the threshold t according to Kapur's method [69], as in Equation (17).The optimal threshold δ is obtained by adjusting the threshold t iteratively to maximize the upper bound of the a posteriori entropy.When the threshold δ is determined, damaged roofs will be identified following Equation (6).

Study Area
The study area (Figure 7) is located in the area surrounding Haiti's National Palace, Port-au-Prince, Haiti.The site is approximately 776,660 square meters of flat terrain with dense buildings.On Tuesday, 12 January 2010, Haiti was hit by an earthquake that struck nearly 15 miles from Haiti's capital, Port-au-Prince.According to official estimates, 316,000 people were killed, 300,000 were injured, 1.3 million were displaced, 97,294 houses were destroyed, and 188,383 houses were damaged in the Port-au-Prince area and in much of southern Haiti [71].Many buildings in the study area were destroyed by the earthquake.

Data Source
The airborne LiDAR data were acquired between 21 January and 27 January 2010, and had an average point cloud density of approximately 3.4 points per square meter [72], as shown in Figure 8a.A 2D GIS vector data of building footprints of the area, which contains 1875 buildings, as shown in Figure 8b, were provided by a third party.

Data Source
The airborne LiDAR data were acquired between 21 January and 27 January 2010, and had an average point cloud density of approximately 3.4 points per square meter [72], as shown in Figure 8a.A 2D GIS vector data of building footprints of the area, which contains 1875 buildings, as shown in Figure 8b, were provided by a third party.

Study Area
The study area (Figure 7) is located in the area surrounding Haiti's National Palace, Port-au-Prince, Haiti.The site is approximately 776,660 square meters of flat terrain with dense buildings.On Tuesday, 12 January 2010, Haiti was hit by an earthquake that struck nearly 15 miles from Haiti's capital, Port-au-Prince.According to official estimates, 316,000 people were killed, 300,000 were injured, 1.3 million were displaced, 97,294 houses were destroyed, and 188,383 houses were damaged in the Port-au-Prince area and in much of southern Haiti [71].Many buildings in the study area were destroyed by the earthquake.

Data Source
The airborne LiDAR data were acquired between 21 January and 27 January 2010, and had an average point cloud density of approximately 3.4 points per square meter [72], as shown in Figure 8a.A 2D GIS vector data of building footprints of the area, which contains 1875 buildings, as shown in Figure 8b, were provided by a third party.

Damaged Roof Detection
In this paper, the detection of damaged roofs is based on damage interpretation using post-event airborne LiDAR data.The process includes three key procedures: data preprocessing, feature extraction and damaged roof detection.Five algorithm parameters, as shown in Table 1, are estimated experimentally or automatically during the workflow, which will be discussed in the Section 4.3.During data preprocessing, 1875 individual buildings are extracted from the raw experimental airborne LiDAR data guided by building footprints, and a dense DSM of each building is generated with a grid cell size of λ from the building points.During feature extraction, first, contour clusters are extracted from each individual building's contours, which are generated based on the contour interval ε from the dense DSM.Two examples of the procedure are shown in Figure 3. Second, the shape chaos index of each contour cluster is calculated using the shape similarity measurement based on the normalized FDs.Contours of the contour cluster are clustered into some groups using the shape similarity threshold ω.Entropy is calculated based on the probability of each contour group and is normalized as the shape chaos index.Third, the 3D shape descriptor of each roof is calculated based on the shape chaos indexes and the area weights of the contour clusters.In damaged roof detection, damaged roofs are discriminated from intact roofs using the optimal threshold δ, which is automatically selected using the maximum entropy threshold based on the roof histogram with a bin width of φ, as shown in Figure 9.The results of damaged roof detection are given as building footprints in a polygon shape file, as shown in Figure 10.

Accuracy Evaluation
For validation, the results are compared to manually labeled reference data based on remote-sensing-based building damage assessment data on the 2010 Haiti earthquake [73], where buildings were classified using the EMS-98 criteria.In these criteria, safe buildings with intact roofs and damaged walls are given grades of 1-3, and heavily or completely collapsed buildings with damaged roofs are given grades of 4-5 [74].To evaluate the accuracy of damaged roof detection in this paper, buildings with grades of 1-3 are re-categorized as intact roofs, and buildings with grades of 4-5 are re-categorized as damaged roofs.
The results of the validation were compared to the reference data, and the accuracy indices, including the overall accuracy (OA), Kappa accuracy (KA), completeness (Compl.)rate and correctness (Corr.)rate, are listed in the confusion matrix in Table 2.The results show that the damaged roof detection technique performed well at classifying both intact and damaged roofs, as shown in Figure 10: the detection technique correctly identified 652 out of 767 (85.01%) damaged roofs and 985 out of 1108 (88.90%) intact roofs.The overall accuracy is 87.31%, and the Kappa accuracy is 73.79%.Some examples of typical buildings, including intact, surface-damaged and structurally damaged roofs, are shown in Figure 11.A surface-damaged roof with debris (Figure 11a) and intact roofs including a flat (Figure 11b), hipped (Figure 11c), gabled (Figure 11d) and pyramid (Figure 11e) roof are correctly identified; this result can also be reliably obtained using other surface damage detection methods using features based on points or planar segments.A completely collapsed building with intact walls as shown in Figure 3e is also detected as a roof with surface damage.Structurally damaged roofs with large planar segments are correctly identified; identifying these structural damages using planar-segment-based features can be difficult.The inclined roof shown in Figure 11f perhaps is a combination of middle and lower story pancake collapse.The inclined roof can be easily confused with an intact single sloping roof; therefore, this structural damage is difficult to identify using segment-based detection methods but is easily identified by the method based on the 3D shape descriptor.Top story pancake collapse, as shown in Figure 11g, with a twisted roof is distinctly represented by contour clusters and is correctly identified.Some intact regular arched roofs are also correctly identified, as shown in Figure 11h, i.
Figure 9 presents the distributions of intact and damaged roofs using a histogram based on the 3D shape descriptor.The undetected damaged roofs mainly show distributions in the critical region below the threshold, whereas undetected intact roofs are widely distributed above the threshold.The main cause of damaged roofs being undetected is that the number and extent of contours affected by damages in some partially damaged roofs are relatively low, which results in the 3D shape descriptors having small values.For example, there is a gap caused by damage on the roof, as shown in Figure 12a.The gap containing points is represented as a contour cluster with various shapes, where the outermost contour has the maximum elevation and is easily detected by the shape chaos index based on contour clusters.However, the area of the gap is very small compared to the whole roof; therefore, the 3D shape descriptor is understated.There are three main reasons for some intact roofs being undetected.The first reason is that the 3D shape characteristics of many shanties exhibit similarities with those of damaged shanties.The second reason is that small components or sundries on intact roofs lead to various irregular contours, as shown in Figure 12b.The third reason is that some intact roofs have regular but dissimilar contours, such as L roofs, as shown in Figure 12c.The shape transformation among these contours is not a translation or scaling transformation; rather, it is an approximate distance transformation.As a result, the normalized FD with invariance under shifting and scaling produces large values of the corresponding shape similarities as well as the 3D shape descriptors.Therefore, the normalized FD is not suitable for measuring the shape similarities among these contours.Some intact roofs are misidentified as damaged because the 3D shape descriptors are overstated for the above reasons.the 3D shape descriptors.Therefore, the normalized FD is not suitable for measuring the shape similarities among these contours.Some intact roofs are misidentified as damaged because the 3D shape descriptors are overstated for the above reasons.In addition, many roofs had already experienced substantial decay prior to the earthquake and were often unfinished [75], thus incorrectly appearing as earthquake-induced damage in post-disaster data.Nevertheless, in the context of earthquake emergency response, undetected damaged roofs are considerably more problematic than misidentified damaged roofs [42].

Parameter Selection and Sensitivity Analysis
In damaged roof detection based on the proposed 3D shape descriptor, four parameters are experimentally estimated, as shown in Table 1.Several key parameters, such as the grid cell size λ, the contour interval ε and the shape similarity threshold ω for contour clustering, are determined using the characteristics of the airborne LiDAR data.It is generally impossible to have each point of an unorganized point cloud associated to one grid height.Consequently, dense DSM is used to approximate the complex building surface represented by the unorganized point cloud.Considering the size of the roof structure and the computational efficiency, the grid cell size of each individual building's DSM is set as 0.1 m in the experiments.The quality of contours derived from a DSM is associated with the grid cell size of the DSM and the average slope in the topography [76].To achieve a comparable quality of contours derived from a DSM, the average horizontal distance between contour lines of the same vertical interval should be approximately equal to the grid cell size.Therefore, we suppose that the optimal average horizontal distance is the grid cell size λ, and the contour interval can be determined as Equation (18).
ε " tanpθq ˆλ (18) where ε is the contour interval, θ is the average roof slope, and g is the grid cell size of the DSM.
According to general building structure knowledge, the average slope is usually less than 45 ˝; therefore, the optimal contour interval should be no more than the grid cell size of the DSM.To select the key parameters ε and ω, a sensitivity analysis is conducted using the variations in each parameter in a reasonable range while the other parameter is fixed at the value that leads to an optimum accuracy.
The overall accuracy (OA), completeness (Compl.)rate and correctness (Corr.)rate are used to evaluate ε and ω.
The quantitative analysis results for the two parameters are shown in Figures 13 and 14 which compare several reasonable values of ε and ω.As the parameter ε increases, the overall accuracy (OA) and the completeness (Compl.)rate increase at first and then decrease; the correctness (Corr.)rate follows the opposite trend.When the parameter ε is 0.08 m, which is less than the DSM's grid cell size of 0.1 m, the OA reaches the local maximum.The result verifies that a relatively reliable result can be achieved when the contour interval is no more than the grid cell size of the DSM.The exact roof slope estimated based on the point cloud can be used to select an ideal contour interval for each building.The OA, Compl.and Corr.also increase at first and then decrease as the parameter ω increases.When the parameter ω is 0.02, the OA reaches a local maximum.It is difficult to quantitatively analyze the parameter ω because it is affected by many factors such as the roof structure, the quality of the contours and the characteristics of the shape similarity measurement based on the normalized FDs.The sample training of contour clusters represents a good option toward selecting a proper parameter ω.Based on the experiments, the values of ε and ω are set as 0.08 m and 0.02, which produce the optimum accuracy in damaged roof detection, as shown in Table 2.

Comparison
The main objective of this paper is for reliable and complete damaged roof detection using post-earthquake airborne LiDAR point clouds.The relevant methods in the literature usually only focused on one or several damage types of roofs.Consequently, in this section, the damage types and the accuracy of the proposed method are compared to those of other methods of damaged roof detection using 3D geometric features extracted from post-earthquake data.
Buildings are usually characterized by various damage types in a serious earthquake zone.In this situation, if a method can detect more damage types of buildings, it will have wide applicability, especially for earthquake emergency response.According to the damage types mentioned in Section 1.1, the damage types of roofs that can be detected by building damage detection methods using 3D geometric features are listed in Table 3.
Table 3.The applicability comparison results and damage types that can be detected by the building damage detection methods using 3D geometric features.

Comparison
The main objective of this paper is for reliable and complete damaged roof detection using post-earthquake airborne LiDAR point clouds.The relevant methods in the literature usually only focused on one or several damage types of roofs.Consequently, in this section, the damage types and the accuracy of the proposed method are compared to those of other methods of damaged roof detection using 3D geometric features extracted from post-earthquake data.
Buildings are usually characterized by various damage types in a serious earthquake zone.In this situation, if a method can detect more damage types of buildings, it will have wide applicability, especially for earthquake emergency response.According to the damage types mentioned in Section 1.1, the damage types of roofs that can be detected by building damage detection methods using 3D geometric features are listed in Table 3.Many methods [40][41][42]44] extract debris for damaged roof detection based on the planarity of the roof.These methods are suitable for detecting surface-damaged roofs but present difficulties in detecting structurally damaged roofs because there is minimal debris in a structurally damaged roof.Vetrivel et al. [29] only focused on the gaps caused by damage and could not detect these building damage types.Gerke and Kerle [32] detected surface-damaged and inclined buildings based on the spatial relation between observable features such as façades, roofs and rubble piles.However, their method cannot detect inclined planes (1) and middle and lower story pancake collapse (4a, 4b, 5a, 5b) because façades are intact and because rubble piles do not exist for these damage types.Shen et al. [31] was able to detect structurally damaged roofs using the geometric axis lines of the roofs.However, it was difficult to detect inclined flat roofs because this building damage type is easily confused with an intact single sloping roof.Fernandez Galarreta et al. [4] was able to detect both surface-and structurally damaged roofs via visual analysis.
A fair comparison between the precision of these methods depends on the type and accuracy of the remote sensing datasets, the number of damage classes, and the type and accuracy of the reference data.The methods in [42,44] and the proposed method detected intact and damaged roofs from different study areas using airborne LiDAR point clouds of the Haiti earthquake, and the results were evaluated using a damage assessment performed by the Global Earth Observation-Catastrophe Assessment Network (GEO-CAN).Labiak et al. [42] obtained an overall accuracy of 73.40% and a Kappa accuracy of 27.51%.Oude Elberink et al. [44] achieved an overall accuracy of 60%.In contrast, the proposed method achieved a high overall accuracy and Kappa accuracy of 87.31% and 73.79%, respectively.
The comparison suggests that the proposed method based on the 3D shape descriptor automatically detects both surface-and structurally damaged roofs and achieves a higher accuracy in damaged roof detection using the 3D geometric features extracted from post-earthquake airborne LiDAR point clouds.This is because the 3D shape descriptor provides a comprehensive description of both surface and structure features based on the shapes and spatial relations of building contours, thus significantly improving the completeness and accuracy of damaged roof detection.A major drawback of the proposed method is that it is difficult to automatically identify specific damage types using the 3D shape descriptor.However, the first priority is to reliably and completely detect damaged roofs in earthquake emergency response.

Conclusions
This paper focuses on both surface-damaged and structurally damaged roof detection using post-event data and proposes a 3D shape descriptor based on building contour clusters derived from airborne LiDAR point clouds for practical application.
The novelty of the 3D shape descriptor is that the 3D shapes of complex roof surfaces are quantitatively described through the spatial patterns of contours.The significant 3D features of surface and structural damages are characterized as a chaotic spatial pattern of contours and are represented as a group of contours with chaotic shapes.The shape chaos index of contour clusters quantitates the chaotic spatial pattern as the information entropy of shape similarities among contours by integrating the spatial attributes and spatial relationships of complex roof surfaces at the regional level.The 3D shape descriptor quantitatively describes the characteristics of chaotic 3D shapes of surface and structural damages at the whole roof level by combining the shape chaos indexes using the area weight.
A significant performance improvement is achieved through the use of the novel 3D shape descriptor based on contour clusters compared to other geometric features directly extracted from point clouds.The experiments produce good results for the airborne LiDAR point cloud used in Haiti earthquake damaged roof detection, which classifies roofs solely using the 3D shape descriptor with the maximum entropy threshold.The overall accuracy and Kappa accuracy are 87.31% and 73.79%, respectively.This damaged roof detection approach is quite valuable for rescue efforts in earthquake-struck areas, especially when pre-event 3D data are difficult to obtain.Future work will be devoted to identifying specific damage types of damaged roofs detected using the proposed 3D feature, as well as to improving the shape analysis between contours to increase the robustness of the 3D shape descriptor.

Figure 2 .
Figure 2. Flow chart of damaged roof detection.

Figure 3 .
Figure 3.The data processing samples: (a) The image of an intact building with an arched roof; (b) The airborne LiDAR point cloud; (c) The dense DSM; (d) Contour clusters colored with different colors; (e) The image of a completely collapsed building with intact walls; (f) The airborne LiDAR point cloud; (g) The dense DSM; (h) Contour clusters colored with different colors.

Figure 4 .
Figure 4. Three building surface regions belonging to different structures are segmented by three contour clusters and are different colors: (a) Contours of the building; (b) Contour clusters of different structures that are represented as nodes colored with the same color in a contour tree.

Figure 5 .
Figure 5.The spatial patterns of complex roof surfaces are represented by contour clusters: (a) The complex surface of Haiti's National Palace, which was destroyed in the 2010 Haiti earthquake (obtained from Google Earth); (b) Contour clusters of the roof are rendered in different colors according to different surface or structural morphologies, such as contours of surfaces with debris (red), twisted surfaces (purple) and inclined roof structures (blue), which are jagged or irregular, and contours of intact surfaces and structures (green) are compact and regular.

Figure 6 .
Figure 6.The complete workflow of the feature extraction algorithm.

Figure 7 .
Figure 7.The location of the study area.

Figure 7 .
Figure 7.The location of the study area.

Figure 7 .
Figure 7.The location of the study area.

Figure 8 .
Figure 8. Data sets of the study area; (a) The airborne LiDAR point cloud; (b) The 2D GIS vector data of building footprints.

Figure 9 .
Figure 9.The distribution of all intact and damaged roofs based on the 3D shape descriptor and the shape similarity threshold δ.The damage information of roofs is from the reference data.

Figure 10 .
Figure 10.The validated results of damaged roof detection, where correctly detected intact roofs (identified intact roofs) are colored green, correctly detected damaged roofs (identified damaged roofs) are colored red, undetected damaged roofs (misidentified intact roofs) are colored yellow, and undetected intact roofs (misidentified damaged roofs) are colored blue.

Figure 12 .
Figure 12.Contour clusters of typical buildings that are difficult to identify using the proposed 3D shape descriptor: (a) Partially damaged roof with gaps; (b) Intact roof exhibiting complex features; (c) Intact L-shaped roof.

Figure 12 .
Figure 12.Contour clusters of typical buildings that are difficult to identify using the proposed 3D shape descriptor: (a) Partially damaged roof with gaps; (b) Intact roof exhibiting complex features; (c) Intact L-shaped roof.

Figure 13 .
Figure 13.Sensitivity test of the contour interval ε.

Figure 13 .
Figure 13.Sensitivity test of the contour interval ε.

Figure 13 .
Sensitivity test of the contour interval ε.

Figure 14 .
Figure 14.Sensitivity test of the shape similarity threshold ω.

Figure 14 .
Figure 14.Sensitivity test of the shape similarity threshold ω.

Table 1 .
The parameter list for damaged roof detection based on the 3D shape descriptor.

Table 2 .
Confusion matrix assessing the accuracy of damaged roof detection based on the 3D shape descriptor.

Table 3 .
The applicability comparison results and damage types that can be detected by the building damage detection methods using 3D geometric features.