Building outline extraction from als point clouds using medial axis transform descriptors

Automatic building extraction and delineation from airborne LiDAR point cloud data of urban environ- ments is still a challenging task due to the variety and complexity at which buildings appear. The Medial Axis Transform (MAT) is able to describe the geometric shape and topology of an object, but has never been applied for building roof outline extraction. It represents the shape of an object by its centerline, or skeleton structure instead of its boundary. Notably, end points of the MAT in principle coincide with corner points of building outlines. However, the MAT is sensitive to small boundary irregularities, which makes shape detection in airborne point clouds challenging. We propose a robust MAT-based method for detecting building corner points, which are then connected to form a building boundary polygon. First, we approximate the 2D MAT of a set of building edge points acquired by the alpha-shape algorithm to derive a so-called building roof skeleton. We then propose a hierarchical corner-aware segmentation to cluster skeleton points based on their properties which are the so-called separation angle, radius of the maximally inscribe circle, and deﬁning edge point indices. From each segment, a corner point is then estimated by extrapolating the position of the zero radius inscribed circle based on the skeleton point positions within the segment. Our experiment uses point cloud datasets of Makassar, Indonesia and EYE-Amsterdam, The Netherlands. The average positional accuracy of the building outline results for Makassar and EYE-Amsterdam is 65 cm and 70 cm, respectively, which meet one-meter base map accuracy criteria. The results imply that skeletonization is a promising tool to extract relevant geometric information on e.g. building outlines even from far from perfect geographical point cloud data. an the CC


Introduction
Mapping building roof outlines, also called building footprints, is essential for digital base map cartography, planning, surveillance, infrastructure management and sustainable city design. Several urban-related applications such as cadaster maintenance and building taxation require building outlines at a routine basis. Extracting building boundary lines manually is expensive and time consuming, especially in urban scenes. Research on extracting building outlines automatically from high-resolution data remains challenging due to the complexity of roof structures and variations in the design of our urban environment. Up to now, building out-lines are typically digitized from multiple aerial images. The use of aerial images is preferred above other data sources as human operators can easily detect buildings on such images. For automation purposes, image disadvantages such as shadows, trees covering building roofs and color variations may increase the extraction error. Moreover, relief displacement may cause problems when using an orthoimage to obtain building boundaries in case of unfavourable image acquisition angles [61] . Airborne Laser Scanner (ALS) point cloud data is an alternative data source. ALS point clouds have been used as a major data source for mapping applications for a few decades [2] . The ability of ALS point clouds to provide many accurate and undistorted 3D points makes it suitable as data source for object extraction. Man-made urban objects (buildings, roads, canals) typically have symmetric shape with straight lines and sharp corners. Such characteristics enable automatic boundary outline extraction from an ALS point cloud. with an efficient algorithm is expected to provide better building outlines.
Medial axis transform (MAT), is a powerful shape extraction technique that provides a compact geometrical representation while preserving topological properties of the input shape [1,4] . The MAT was introduced by Blum [5] to describe biological shapes. Since then, it has been used for applications in image processing and computer vision. However, MAT has a fundamental drawback, which is its instability to small perturbations of the input shape, which then may disturb the topology of the MAT branches [3,9] . Moreover, wider deployment of MAT to extract shapes analysis from surveying quality data with its associated problems, is still challenging [4] .
In principle, the MAT can be implemented for urban mapping purposes, particularly to extract building shapes by detecting its corners. As illustrated in Fig. 1 , corner points (red) are detected when the tip of a skeleton branch (blue line) touches the building boundaries (black). However, generating a MAT skeleton from point clouds is a challenging problem as such data contain fuzzy borders, in particular, when data is missing [25] . This makes the application of MAT for airborne point clouds difficult.
Corners are important local features and knowledge on their location can minimize further data processing without losing specific features of the original object shape [35,36] . Given an airborne point cloud of an urban area, we propose a method for extracting building outlines automatically by detecting accurate roof corner points based on MAT descriptors.

Related work
Work on the development of building outline extraction from various remote sensing data has been intensified in parallel with the increased interest in GIS (Geographic Information System) digital map products [49] . Combining different data sources to extract building outlines is believed to increase the detection rate and accuracy compared to using a single data source as more features can be used. Nevertheless, fusion of data of different sensor types is not a trivial task as fusion is hampered by dissimilar resolution, alignment issues or mismatches in feature information caused by sensor characteristics or differences in viewpoint during acquisition [60] . Manno-Kovacs and Sziranyi [54] proposed an Orientation Selective Building Detection method to detect buildings from a combination of aerial and high-resolution satellite images. They apply active contouring for obtaining smooth and accurate building outlines. However, inhomogeneous buildings are sometimes only partially detected and any object connected to the building (e.g. trees) can result in false positive detections. Zhao et al. [55] , Awrangjeb [56] , and Li et al. [57] combined LiDAR point clouds and aerial images to detect buildings and obtain smooth building outlines by regularization and mathematical morphology. Building outline errors occured due to failures to determine the building principal directions during regularization [55,56] or line redundancies after simplification [57] especially, in case of low point density. Our previous study [37] proposed an extended Hough transform method using ordered lists of points to detect building boundary segments from airborne point cloud data. Hierarchical filtering is applied to remove spurious lines. That method outperformed existing state of the art methods in terms of correctness and quality metrics on benchmark dataset. However, the method is likely to introduce false corners, especially for buildings of complex shape, as spurious lines may still exist in the final step.
Various definitions of MAT or skeleton found in literature correspond to different methods for computing the MAT leading to different results with different properties. In general, MAT algorithms typically focus on deriving the geometric location of the centerline or medial axis of a surface, so-called skeletonization [11] . Up to now, numerous skeletonization methods and their application for 2D and 3D object description are available in literature. The existing skeletonization methods are often categorized into four main approaches: -morphological thinning-based methods that was first applied for discrete binary images [17] , and improved by Huang et al. [22] and other related notions on the development of 3D thinning algorithm [23,48] . -geometry-based methods using medial axis transformation for planar shape [18] including Voronoi diagrams [19] and Delaunay triangulation [47] . -distance-based functions such as skeleton generation and centerline by the distance transform [40] , skeletonization by discrete Euclidean distance maps [41] , and Euclidean skeleton based on connectivity criterion [42] . -general-field functions which are generated by functions rather than use distance function, for example by replacing the nonlinear distance with a linear transform [16] , using Newtonian potential model to replace the distance function [43] , and using Electrostatic Field Theory (EFT) function to generate potential distribution inside the object [44] .
Reviews on skeletonization methods and its applications have been discussed by Saha et al. [12,13] , Pavlidis [21] , and Amenta et al. [24] . As many skeletonization algorithms were designed more for image analysis, we limit the scope of our study on skeletonization for point cloud data.
Several works on MAT using geographical data for various purposes have been conducted. Haunert and Sester [30] applied straight skeleton extraction to derive linear representations of polygons and road centerlines from a cadastral dataset. Yirci et al. [29] extracted detailed pedestrian networks by generating a centerline using two skeleton operators (straight skeleton by parallel thinning and medial axis by Voronoi diagram). Methods for river centerline extraction based on Delaunay triangulations remain challenging for certain complex situations e.g. a scenario with a skeleton branching in different directions [33,34] . Broersen et al. [26] used the 2D skeleton of a Voronoi diagram and the 3D skeleton of a shrinking ball for identifying watercourses and deriving its centerlines from classified aerial point clouds. Widyaningrum and Lindenbergh [38] extract the road network of an urban area from a colored point cloud using parallel thinning skeletonization [39] . However, a further generalization step maintaining the road topological order is required for smoothing jaggy skeletal lines yielded by the parallel thinning algorithm.
Ma et al. [6] estimate a 3D medial axis point using a shrinking ball approach based on nearest neighbors and normals from a given set of points. Their work is claimed as a faster and simple method that approximates the 3D MAT based on maximally inscribed balls tangent to two surface points whose center is positioned on a normal line. The ball centers are the medial axis points. The shrinking ball algorithm is not only accurate and computationally efficient but is also considered as the most simple and fast existing surface-skeletonization method [4] . This method is considered as suitable for the geographical case as it is point based, simple, fast, and scalable [3] . However, to use this method, fine sampling is required to directly obtain a high quality skeleton. The fact that the shrinking ball approach can only result in unstructured skeleton points while disregarding the topology of the skeleton branches makes this algorithm not directly applicable for applications involving surveying data. For use in practical applications, the MAT consisting of medial points and corresponding maximally inscribed circles (in 2D) or balls (in 3D) needs further processing [10,14,15] .

Methodology
Our research focuses on the adaptation of MAT for extracting building outlines from noisy point cloud data required for mapping and spatial modeling purposes. We extend the work on the iterative shrinking ball algorithm and develop a strategy to exploit skeleton features to accomplish the goal of accurate building outline extraction. A new approach for skeletal point segmentation is also proposed in this study. The proposed method achieves stateof-the-art in handling noisy surface boundaries and reconstructing the building outlines. It requires minimal human interaction by optimizing the use of skeleton-based features. Specifically, the contributions of our work are as follows: Overall, our method overcomes some traditional pitfalls of using MAT techniques in case of noisy input.
As this study requires using the MAT in 2D space, we adapt the 2D shrinking circle algorithm by Ma et al. [5] . The general workflow of our proposed method for automatically extracting building roof outlines consists of four main steps (see Fig. 2 ). First, building boundary points are extracted by an alpha-shape algorithm [45] . Next, the boundary points are transformed into its 2D MAT or skeleton points using the 2D shrinking circle algorithm. Third, we then apply our MAT segmentation to segment the MAT points by exploiting their geometric attributes. The segments are then used to detect corner points. Fourth, polygonization is carried out to form a 2D closed polyline based on the detected corner points.
This study uses the extended shrinking circle approach that implements the denoising heuristic as proposed by Peters [3] . We define the skeleton of an object surface S as a set of center points c of maximally inscribed circles B ( c, ρ) in S (see Fig. 3 ) where ρ denotes the radius of such circle.
The 2D skeleton points are also called medial axis points. By associating the circle radius ρ function to the set of medial axis points, we obtain the so-called Medial Axis Transform (MAT). As shown in Fig. 3 , the medial axis points (red points) form the MAT skeleton (blue lines) of a rectangular object S . Each maximally inscribed circle (in grey) touches at least two points of the boundary of S (black outline). Center points of any circle that is not maximal or not inscribed in S (green circles) are dismissed and not considered as medial axis points.
To provide a clear and coherent narrative, further details on each methodological step are provided in the following: Section III.1 describes the alpha-shape algorithm. Section III. 2 and III. 3 provide necessary details on the shrinking circle method and skeletal points extraction, respectively. Skeletal point segmentation is described in Section III. 4. The last step, corner point estimation and building outline generation, is explained in Section III. 5. Evaluation methods for building outline extraction used in this research are discussed in Section III. 6.

Alpha-shape
Given the segmented building points, creating building outlines starts with boundary point selection by the alpha-shape algorithm as introduced by Edelsbrunner [45] . An alpha-shape is well known for its capability to preserve small shape details of a finite point set at a required level of detail. The 2D alpha-shape is constructed based on the 2D Delaunay triangulation of the input points. The method identifies boundary points that are defined in terms of a parameter α ≥ 0, which controls the level of detail of the boundary shape. Given a set S of points on a plane and a value of α, the algorithm works as follows: 1. Compute the Delaunay triangulation DT( S ) of S. All edges in DT( S ) are candidate for the alpha-shape S α . 2. For all edges e of DT( S ) with end points p and q , say: a. Find two circles B pq 1 and B pq 2 of radius α with center c pq (1) and c pq (2) containing end points p and point q of the same edge e . The circles are defined in terms of the below circle centers: Where e is the length of the edge between end points p and q. b. If at least one of the circles contains no points from S in its interior, e is a valid boundary edge ( Fig. 4. a), otherwise the edge is removed ( Fig. 4. b). 3. The union of all valid boundary edges forms the alpha-shape S α ( Fig. 4. c) The value of α is a real number with 0 ≤ α ≤ ∞ . As α approaches 0, the shape may shrink, develop holes and may become disconnected. In the extreme case, the value of α = 0 results in the data points itself. When α increases towards infinity, the alpha shape approaches the convex hull of the set S of points. In case of geospatial point clouds, the point density is often varying, and is depending on sensor characteristics and measurement geometry and an appropriate value of α should be chosen accordingly.
To identify boundary points of unordered building roof points in our study areas, we decided empirically for an α-value between 0.3 and 0.5

The shrinking circle principles
Given are a set of noisy edge points V on a surface S with corresponding normal vectors N . The MAT points are defined as the set of centers c and corresponding radius ρ of maximally inscribed circles B ( c, ρ) in S that are bi-tangent to the boundary S . The circle B and corresponding circle center c are denoted as medial circle and medial axis point, respectively.
The basic principles of the shrinking circle method (see Fig. 5 ) are as follows: 1. A medial circle touches the surface in at least two points ( p, q ) where p , q ∈ S . 2. Following the line defined by normal vector N p of edge point p , the radius ρ of a circle B p decreases iteratively until B p touches S at q , where q = p and the circle center c is on the line through N p . Iteration stops if the maximal B p circle is found. 3. A medial circle is a maximal empty circle, which means it contains no surface points.

Skeletal points extraction
To obtain the MAT of surface S , medial axis points c ( p ) are computed. Hence, the maximal inscribed medial circle B for all sample points p in S is computed by the following steps: 1. An initial circle B init of p is defined based on an initial radius ρ init . The ρ init value is set sufficiently large e.g. equal to the largest distance between two input points.
2. Given ρ k p , where k = {1, 2, .., i } denotes the k-th iteration step, the circle center c k p is given by: 3. Find the surface point q k p ∈ S closest to c k p such that q k p = p . 4. Test for circle maximality for the circle defined by q k p and p : a. If the distance from c k p to q k p equals the radius of the circle ρ k p , the circle B k p is maximal and c k p is a medial axis point. b. Otherwise, compute the radius of the next shrinked circle ρ k +1 p using the following equations: The iteration will stop when the medial axis point as described in step 5.a. is found. Fig. 5. a. shows consecutive shrinking of a circle touching S at point p, which results in a medial circle B i p and a medial axis point c i p in the last iteration. Given a defined inside and outside of surface S , the MAT consists of two components: one part inside surface S (N p ) , consisting  of the so-called inner medial axis points, and another one outside surface S (-N p ) , the outer medial axis points.
For each p ∈ S , the corresponding inner and outer MAT points are computed by iterating steps 2 to 4. The inward normal N p is used for the inner MAT calculation, while the outward normal-N p is used for the outer MAT calculation. Fig. 5. b. shows the geometry for calculating the medial axis point c p and the direction of the normal vector N p for the inner circle (black arrow) and outer circle (red arrow).
Noise handling is an essential step to overcome the sensitivity of MAT to noisy boundaries. In case a small bump or noise ex-ist on the input surface, a circle may get shrunk too much which then likely results in undesireable medial axis points. Such overly shrinked circle, typically has a small separation angle α. The separation angle α (see Fig. 5 The denoising heuristics technique presented by Peters [2018] is used to select the so-called good circle that is often computed during the above described shrinking approach. The good circle is defined as the last circle in the shrinking approach with a separation angle α k bigger than the separation angle threshold α min .
After this step, every medial axis point c p comes with a number of attributes. Attributes of each MAT point m p are medial axis point c p coordinate, radius ρ, separation angle α, indices of surface point p and q , and normal vector N p or -N p . Theoretically, using these MAT attributes, the geometry of S can be reconstructed completely.

MAT point segmentation
MAT attributes provide rich information that can be used to group MAT points into different medial segments or branches. In our case, segmenting the points into different branches is used for detecting the corner points. One issue when obtaining a skeleton from a point cloud with the shrinking circle algorithm is that the point cloud provides an unstructured point sampling of S . Also for the MAT points resulting from steps 1 to 4, adjacency relations are initially not known. For further application of the MAT, two useful observations to identify MAT point connectivity are as follows: -MAT points heading towards the same turning point or corner are considered as one segment. Fine sampled surface points of a square shape S in Fig. 6. a result in fine MAT points in which some of the MAT points gradually approach a specific turning point. In this sense, each MAT point created from a maximally inscribed circle, touching at surface point p and q appoint to a turning point that is equally located between surface point p and q . As illustrated in Fig. 6. b, the median value of two surface points p and q (in red text) is similar to the corner's index (76).
-MAT points are expected to have a separation angle α close to 90 °As shown in Fig. 6. c, MAT points of a rectangular shape have separation angles that are distributed around 90 °I n practice, surface points are not perfectly distributed and noisefree and are not as regular as shown in Fig. 6. a. Small perturbations on the surface boundary create so-called skeletal noise [7,8] . When detecting shape corners, skeletal noise may induce false segments, which then results in false corners. Our segmentation criterion relies on the proximity of the turning point location of the boundary point. Intuitively, MAT points that lie close to each other and have similar features values are grouped together. Moreover, we expect that the radius ρ will gradually change along a segment branch.
Based on aforementioned observations, we use three global thresholds and four MAT-derived features for segmenting the MAT points. The global thresholds are not related to MAT and are defined to increase the segmentation accuracy. The global thresholds are: G.1. A buffer distance bf from the object surface S . Only a MAT point located within the specified buffer will be considered for segmentation. This threshold is used to exclude unusable outer MAT points resulting from the maximal circle of two edge points with outward normal. These MAT points are typical noise located far from the surface points. G.2. Minimum number of points for each segment minPts. Any segment having less points than the given minPts is considered as segment noise .

G.3. Point index interval
pt sets the minimum distance between two candidate corner points, as expressed below: In Eq. (7) , l is the minimum required edge length and r is the point cloud interval. The point index interval pt criterion is designed to avoid having false or extra corners at a certain minimum determined edge length in case of short and noisy boundaries. For example, given a set of points with 0.5 m point cloud interval r , we require to extract building edges of minimum length l = 2.5 m, thus, pt is set to 3. Imagine that point 13 in Fig. 7 has the same medial properties as point 11, 16, and 20. Point 13 will not be considered as a corner point as it has less than 3 point difference to point 11.
The customized features derived from the MAT attributes, or MAT-derived features are described as follows: The normal angle differences δN between 2 consecutive edge points is used to determine candidate corners K p for p ∈ S . Detection of candidate corners K p is used to obtain first estimate of the number of building corners which is later used to remove false segments in case of noisy edges. At noisy edges, segments with small median difference may be formed which later may result in false corners. Consider two adjacent normal vectors N p i and N p i −1 of point p and A surface point p initiates a candidate corner point K p if the angular differences δN to its two adjacent points ( p i −1 and p i +1 ) are above the given angle threshold δN . That is expressed in Eq. (11 ). 6. If | me d pq − r idx | ≤ pt, then MAT point m p is assigned to medial segment Mseg ( r ). 7. Any medial segment Mseg having less member points than minPts, as defined in G.2, will be removed. This step will eliminate false segments that may be formed in case of flaws on the edges.
Medial segments are used next to estimate real corners where one medial segment corresponds to one corner.

Corner point estimation
Instead of appointing edge points as corners, we rather estimate the position of corners based on the medial axis point positions and their corresponding radius. The radius ρ of the maximally inscribed circles of MAT points will gradually decrease towards corners. Each medial segment ideally contains a set of MAT points with gradually decreasing radii. The location where the radius will become zero ( c ρ=0 ) , typically identifies the location of the corner point corresponding to a medial segment is estimated as follows. Fig. 9 shows how the radius ρ of the MAT point depends linearly on the x coordinate. Therefore, a line is fitted by PCA (Principle Component Analysis) through the ( x, ρ) points of the segment at hand. The x -coordinate corresponding to zero radius ( ρ = 0) of the fitted line L r (as indicated by the blue line in Fig. 9 ) is reported as the x coordinate of the corner point. The ycoordinate of the corner point is obtained in a similar way.
In case of an edge with heavy noise, false segment may remain. A spatial filtering step is necessary as final filter to remove any spurious estimated corner c . This spatial filtering step preserves any corner point that is within a specified radius from the surface points p ∈ S . In our case, we use 1 m as we require building outline result to have positional accuracy at least 1 m. As final step, a closed building polygon is obtained by connecting all corner points consecutively referring to the point indices.

Building outline evaluation metrics
Two different evaluation metrics are applied for evaluating the performance of the proposed workflow in fulfilling the required building outline specification: corner geometric accuracy and corner detection accuracy. To evaluate the geometric accuracy, we compared the position of corners coordinates of the building outline results to the reference [46] . Positional accuracy, also known as geometric position accuracy or location accuracy, is used as main indicator to measure how well the building polygons are positioned with respect to its true position within an absolute georeferenced system. We use the RMSE (Root Mean Square Error) to measure the average of the squared differences between building corner positions (X and Y coordinate) in the reference and in the result. The RMSE of a complete building is calculated for all detected building corners with respect to the position of corresponding reference corners.

RMSEx
Where: X res , Y res = Coordinates of resulting corner points X ref , Y ref = Coordinates of corner points in the ground truth n = total number of corner points Due to the complexity of some building, not all corners may be detected completely. Therefore, we evaluate the corner detection accuracy by means of three retrieval measurements: recall, precision and F1-score [59] . Precision is used to measure the exactness or fidelity, whereas recall is used to measure the completeness. The F1-score is the weighted mean of precision and recall.

P recision =
T P T P + F P (15) Recall = T P T P + F N (16) For this purpose, a corner point is considered a True Positive (TP) if it is located within 1-meter radius from the corresponding corner reference, while any undetected corner including corners with an offset of more than 1 m from the corresponding reference In the building polygon in Fig. 10 ., the number of correct corners (TP) is 4, while the number of false corners FP (inside the green ellipse) is 1, and the number FN of undetected corners or corners with an offset of more than one meter (as indicated by the blue circles) is 2. This configuration gives Precision = 0.8, Recall = 0.67, and its F1-score is 0.73. Fig. 11 summarizes the proposed method to detect building corners from a given cluster of building points ( Fig. 11. d). The parameter thresholds, as discussed above, are set empirically depending on the point density and the required specification. Additional pre-processing is necessary in case clustered building points are not available. In this case, an initial classification and/or semantic segmentation processing step is required. We use a classified ALS point cloud whose points are labelled according to their object class (building, ground, and unclassified). To group points belonging to one building, points are clustered by applying the DBSCAN algorithm [50] . As a result, different buildings will have different cluster number and theirs points are labelled according to the corresponding cluster number. Once the clustered building points (as presented in Fig. 11. d) are available, boundary points need to be extracted ( Fig. 11. e), for which we use the alpha-shape algorithm. The resulting building boundary points are then used as input for the MAT shrinking circle algorithm (blue points Fig. 11. f). MAT points filtering and segmentation is applied based on their separation angle and median index value ( Fig. 11. g). Each medial segment generates a different corner candidate. Positions of the corners are extrapolated linearly using PCA ( Fig. 11. h).

Experiments of the study areas
For the experiment, we use three study areas with different landscape characteristics and airborne LiDAR point cloud specifications. The first dataset represents a sub-urban area of the city of Makassar, Indonesia. The point cloud data was captured in 2012 by a Leica ALS70 instrument and has 7-11 ppm (point per meter) point density. The second dataset is a Dutch national AHN3 point cloud [31] sampling the area of EYE-Amsterdam, the Netherlands. The AHN3 data has a point density of at least 10 ppm and was acquired in 2014. Most buildings in this area have a public or business function. This test set is selected as many of the buildings in this area are considered to have a high complexity in terms of shape and size.
Building points of the Makassar dataset, as presented in Fig. 12. a. in orange, were classified using LAStools [27] . For EYE-Amsterdam, we used the provided building classification of the AHN3 dataset (shown as orange points in Fig. 12 ). The alpha-shape algorithm [28] is then applied to derive the outline of each building. We do not discuss the details of the pre-processing steps further as it combines well-known methods in the field of GIS and remote sensing.
For the Makassar area, the topographic base map scale 1:10.0 0 0 is used as ground truth data. The topographic base map is gen-erated from manual 3D delineation from stereo-images with the same acquisition time as the Makassar airborne point cloud data. For validating the Eye-Amsterdam results, the Dutch building registration dataset BAG ( Basisregistratie Adressen en Gebouwen ) of the 2019 building dataset is used [32] . However, we noticed that several BAG buildings have different shape and size compared to the AHN3 buildings due to different data acquisition time. Thus, the RMSE of the corners is calculated for unchanged building.
In this study, the required positional accuracy for the outline result is at least 1-meter. For the Makassar dataset, the average RMSE for 36 buildings in the study area is ̴ 65 cm, which meets our requirements. There is an exception for one incomplete building (indicated by a red circle in Fig. 11. b.), that has a RMSE of more than one meter. In this case, the building roof is partially covered by dense trees resulting in poor building segmentation. Based on the number of corners from the reference data, the precision, re- Building outline results for the EYE-Amsterdam dataset, as shown in Fig. 12. d, are also good at 0.7 m RMSE on average. For this dataset, there are no classification issues as we use available building points provided by the AHN3. In this case, one of the biggest factors that influences the accuracy of the result is the definition of building roof, in case an overhanging roof exist. This means that the overhanging roof is likely included in a building in AHN-3 data but not in the BAG data, which results in discrepancies. Fig. 13. b shows overhanging roofs on one building in the EYE-Amsterdam area (marked as A in Fig. 13. a). Another issue is found on a building with a curved building outline. The algorithm could not detect points on a curved line, as the normal angle differences δN between edge points is small (less than five degrees). Decreasing the normal angle differences δN threshold would not always mitigate this problem, as it will increase the number of false corners due to noise of the edge points.

General overview
Compared to our previous study on building outline extraction using ordered points aided Hough transform (OHT) [37] , the MAT approach has higher sensitivity to noisy edge points and varia-tions in point density. Small bumps on the edges affect the normal direction of the corresponding point, which later affects the corner detection result. However, the MAT approach has a similar accuracy as our previous work at 0.7 m RMSE for the MAT-based method and 0.57 for the OHT method. Our proposed method is able to handle one of the shrinking ball algorithm limitations that requires the surface normal for each sample point [6] .
The use of the alpha-shape algorithm makes it possible to orient the normals of each edge point automatically when performing inner and outer MAT computation. Inner and outer MAT points together effectively detect all building corners. Fig. 14 shows the ability of our method to detect corners of various building shapes with different point density. Fig. 14 row 1 demonstrates how our algorithm works on sparse building edge points that later result in false segments. The red circle in Fig. 14. b row 1 shows three medial segments around a corner. Two of them are considered as false segments. However, the algorithm delivers the correct corner points only (inside the red circle in Fig. 14. c row 1) even though false segments exist. The method also works in case small perturbations exist on the boundary e.g. due to trees (see orange circles in Fig. 14. a row 3). There is one MAT point (the red point inside the yellow ellipse in Fig. 14. a row 3) resulting from a noisy edge part, that is correctly ignored by the medial segmentation procedure. On the other hand, the algorithm may fail  to extract non-linear shapes like rounded edges (see blue circles in Fig. 14 . row 2) as the algorithm cannot detect significant normal change for boundary points on the rounded edge. In this case, the algorithm can only detect two candidate corners from two medial segments, positioned at both ends of the rounded edge.
As shown in Fig. 15. a, our method is able to obtain building corners (yellow points) that are close to the reference building polygon (green) even when there is a perturbation in the boundary (inside the yellow ellipse). The proposed method, particularly, improves the alpha-shape outline (red). Fig. 15. b compares between corners from our method, the alpha-shape outline result, and the building outline reference.
In addition, we found out that several MAT properties have not been investigated yet, particularly for the outline extraction, e.g. the curvature of consecutive MAT points indicating an asymmetric shape of two edges of different length, may be used for locating but also characterizing corners accurately.

Comparison analysis
In this section, we compare results of our proposed method to those of existing methods on building outline extraction. For   [53] . All these methods use an alpha-shape algorithm to select boundary points, which later result in building corners. The RanSAC-based segmentation and regularization method requires primary building orientations as it regularizes all boundary lines with respect to these orientations and its perpendicular orientations. Another regularization approach, OHT [37] , applies the so-called extended Hough transform on a list of ordered boundary points that enable to detect arbitrary building directions and extract different boundary segments. The comparison metrics considered are the building corners geometric accuracy (RMSE), the computation time of building outline extraction after the boundary points are selected, and the corner detection accuracy in terms of recall, precision, and F1-score.
As shown in Table 1 ., our proposed method has the highest geometric accuracy as well as F-1 score. The RanSAC-based method has higher Recall value as our proposed method, but is likely to have more undetected corners. However, the proposed method has the lowest number of false positive corners. The average computation time for the three methods is considered comparable, although in fact the proposed method is the fastest. However, each method may have different strengths and weaknesses when it is applied on a complex building roof shape. Using the building shown in Fig. 16 . as one example, a reference corner inside the brown circle is detected by OHT, but not by the proposed method. The RanSAC-based method results in two corners with an offset close to 1-meter from the reference. The proposed method fails to detect one building corner (inside the brown circle) due to a wide angle between two consecutive building outline segments that is close to 180 °In this case, the shrinking circle failed to produce a separation angle α close to 90 °needed for corner detection.
Another difference in the results in Fig. 16 . is indicated by the blue circles. The OHT method fails to detect two corners of a short building edge of one meter length, while the RanSACbased method as well as our proposed method successfully detect both corners. All methods produce false corners but the proposed method has the smallest number of false corners.

Computational and complexity analysis
We implemented our method in Python2.7 on an Intel Core 2Duo CPU with 2.4 GHz processors. The computation time for our corner detector algorithm varies from 0.28 to 0.99 s per building, depending on the building size, shape and point density.
Recall that the proposed method has four main steps: boundary point selection by an alpha-shape algorithm, MAT extraction by a shrinking circle approach, MAT points segmentation, and estimation of corners by PCA. However, in case building points are not yet available, an additional classification and segmentation step is required. One of the most common method for 3D point cloud segmentation is seeded region growing. This type of algorithms selects a seed point and adds a point from the neighborhood if it meets a certain criterion. As reported by Shih and Cheng [51] , and Deschaud and Goulette. [52] , the computational complexity of seeded-based region growing for partitioning a point cloud consisting of n 3D points into N segments is O ( n log n ). Rabbani et al. [55] improved plane detection using a smoothness constraint based on the normal angle difference between neighboring points and reported a time complexity of O ( n + N log n ) .
The proposed algorithm works on individual building segments. Suppose a segment contains u number of building points, the computational complexity of our four main steps are as follows: 1. The alpha-shape algorithm for selecting boundary points using Delaunay triangulation, as reported by [58] , has a time complexity of O ( u log u ). If all points of a building segment were boundary points, the computational complexity would be O ( v 2 ) time. This computational complexity for the total procedure of medial axis segmentation is not a problem because the algorithm works just on the boundary points of a single building. Point cloud segmentation is in practice the heaviest computation as it involves the complete 3D point cloud data. For example, in our Makassar test area, the building segmentation step considered all 464.191 points, while the medial axis segmentation considers 43 individual buildings. The number of points per-building in the Makassar test area varies from 96 to 3185 points, that later results in a number of building boundary points varying from 23 to 156 points. Note that the processing of individual buildings can also be parallelized easily.

Conclusions and recommendations
In this study, we have presented a procedure for automatically extracting building outlines from airborne point clouds based on the MAT descriptors generated by the 2D shrinking circle method. Our approach takes advantage of MAT invertibility with its medial axis and the corresponding radius function that allows reconstructing the exact object shape. Building classification is conducted first. A set of building boundary points is then extracted using an alpha-shape algorithm. After applying a shrinking circle algorithm to the input boundary points, MAT points are obtained. To identify corners, we introduce a corner-aware segmentation to group MAT points to their corresponding medial branch. The segmentation combines both global thresholds and several MAT-derived features. Next, the algorithm fits a line to all MAT points of a segment. Based on the corresponding radii of the MAT points, the corner point location is estimated by extrapolating the position where the radius is zero on the fitted line.
The positional accuracy results of the estimated corner points indicate that our method provides a completely new and promising tool for reconstructing the geometric shape of building roofs from scattered airborne point clouds by using the MAT approach. The proposed method performance is highlighted on a number of complex building shapes in and around the EYE building in Amsterdam, The Netherlands. The ability of the proposed method to obtain accurate corners and complete shapes indicates the robustness of our method to small perturbations on the building edge. In case of sparse point intervals, densification of edge points may help to increase the accuracy of the MAT result. Meanwhile, our method has limitations to obtain outlines of an object with a curved or circular shape. In comparison to state-of-the-art methods on building roof outline extraction, our proposed method shows a promising result in acquiring accurate building corners geometrically. Compare to RanSAC and Hough transform based methods, our method is a primitives-free approach that does not require orientation initialization.
Though current skeletonization methods show a progressive development, deployment for wider applications is still challenging. Different applications may require different skeletonization methods and/or MAT descriptors. For future work, we will consider extending the MAT technique for reconstructing other digital map objects, such as road networks, ridges, or streamlines. The application of the proposed workflow for extracting curved lines and for a larger area is also interesting to be investigated further.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.