# Determination of Building Model Key Points Using Multidirectional Shaded Relief Images Generated from Airborne LiDAR Data

Authors: Dong-Cheon Lee; David H. Lee; Dae Geon Lee
Journal: Journal of Sensors
 (2019)

---

## Abstract

Light detection and ranging (LiDAR) data collected from airborne laser scanner system is one of the major sources to reconstruct Earth’s surface features. This paper presents a method for detecting model key points (MKPs) of the buildings using LiDAR point clouds. The proposed approach utilizes shaded relief images (SRIs) derived from the LiDAR data. The SRIs based on the concept of the shape from shading could provide unique information about individual surface patches of the building roofs. The main advantage of the proposed approach is to detect directly MKPs, which are primitives for 3D building modeling, without segmenting point clouds. Depending on the location of the light source, the SRIs are created differently. Therefore, integration of the multidirectional SRIs created from different locations of the light source could provide more reliable results. In addition, the vertical exaggeration (i.e., scalingZ-coordinates) is also beneficial because constituent surface patches of the roofs in the SRIs created with vertically exaggerated LiDAR data are more distinguishable. To determine the MKPs of the roofs, building data was separated from other objects using modified marker-controlled watershed algorithm in accordance with criteria to specify buildings such as area, height, and standard deviation. This process could remove the unnecessary objects such as trees, vegetation, and cars. The curvature scale space (CSS) corner detector was used to determine MKP since this method is robust to geometric changes such as rotation, translation, and scale. The proposed method was applied to simulated and real LiDAR datasets with various roof types. The experimental results show that the proposed method is effective in determining MKPs of various roof types with high level of detail (LoD).


---

## Body

## 1. Introduction

Since the airborne laser scanner (ALS) systems have been commercialized in the late 1990s, light detection and ranging (LiDAR) data collected from ALS systems have been widely adopted as a major source in geospatial information engineering such as city planning, mobile navigation, forest mapping, disaster response, and damage assessment. Because of precise and accurate data acquisition capability, LiDAR is one of the most preferred remote sensing technologies. The advantage of the LiDAR is direct acquisition of dense point clouds with geocoded 3D coordinates over extended areas cost effectively and quickly [1]. LiDAR sensors are not affected by geometric distortions unlike optical sensors (e.g., lens distortion, perspective distortion, and relief displacement). In addition, shadow and variation of the brightness and contrast due to illumination condition affect quality of the optical image while LiDAR data is less influenced by such external factors.Numerous 3D building modeling (or reconstruction) methods have been proposed in the past decades. Existing approaches are classified as data-driven (or bottom-up) and model-driven (or top-down) approaches. The data-driven approach might consider faithful details pertinent to the resolution of the input data, while the model-driven approach could rapidly create building models with appealing appearance [2]. In the model-driven case, building models are recognized by fitting the LiDAR measurements with predefined models. Therefore, a set of predefined models is required as parameterized shapes. If the model library does not include all types of the buildings, reconstruction of the building models might be not completed. In the data-driven case, there is no need to have a prior knowledge of building type. Since parameterization of the building models in model-driven approach is often difficult, the data-driven approach is preferable. The major issue of the data-driven approach is to find the constituent planes and primitives (e.g., edges, vertices, and corners) of the building roofs.Majority of this approach aims to identify and segment individual roof patches that forms the roof shape of a building. Once roof surfaces are segmented, each segmented surface patch is fitted with an appropriate mathematical function. Then, form-lines (or structure lines) and key points are determined by solving the surface equations. Such approach requires topological information, for example, adjacency relationship among segmented surfaces. Therefore, boundaries of the roofs, form-lines of the various roof types, and key points (e.g., corners, intersections, apexes, and vertices) are the most important features for reconstructing 3D building models. 3D building modeling is one of the most prominent applications of LiDAR data [2–5]. Even though there is not a standard procedure of the building modeling available, the common tasks involved with LiDAR data processing include strip adjustment, noise removal or reduction, separation of target objects (e.g., building data extraction) from other features, interpolation (if necessary), and segmentation for 3D modeling. In particular, surface patch segmentation of the roofs is one of the major parts in the building modeling.As one of the earliest researches related to building modeling using LiDAR data, Csathό et al. [6] performed segmentation of LiDAR points by iterative robust sequential estimator to parametrize and organize surfaces. The procedure consists of interpolation, edge detection, seed selection, determination of the best fitting mathematical function, isolate point removing, and gap filling. Rottensteiner and Briese [7] presented a method for building extraction and automatic 3D building modeling. The procedure starts with separation of points on the buildings from those on the other features including terrain using hierarchical robust interpolation of LiDAR points. The initial building regions were determined based on the height difference between terrain surface and digital surface model (DSM). The geometric 3D building models were generated by grouping segmented planar surfaces. Sampath and Shan [8] introduced roof segmentation approach using clustering based on the eigenvalues of the covariance matrix in a small neighborhood and k-mean algorithm. Park et al. [9] proposed surface segmentation method by multilevel and multidirectional 3D chain code that could determine slope and orientation of the individual surface patches. Chen et al. [10] presented adaptive random sample consensus (RANSAC) for roof patch segmentation with building points extracted by morphological filtering and region growing algorithm.Wang and Shan [11] demonstrated building extraction scheme by segmentation of point clouds using local similarity measures based on the discrete computation geometry and machine learning. Yan et al. [12] applied snake algorithm to process LiDAR data that is extended dynamic programming method to address the snake problems with a 2D planar topology using a graph reduction technique to enforce geometric constrains. They could construct complicated types of the roof with topology adjustment and refinement strategies. Sampath and Shan [13] proposed a potential-based approach to estimate number of the clusters while considering both geometry and topology for the cluster similarity. The final step of segmentation separates the parallel and coplanar segments based on their distances and connectivity, respectively. Building reconstruction starts with forming an adjacency matrix that represents the connectivity of the segmented planar segments. Tseng and Hung [14] applied split-and-merge based on octree structure with plane fitting to segment LiDAR points into coplanar clusters.Segmentation is the central issue in most building modeling approaches that use the LiDAR data. Segmentation is to group the point clouds with similar geometric characteristics that could be constituent surfaces of the roofs. Peternell and Steiner [15] separated the LiDAR data into a number of grids. For each grid and its eight neighbors, the surface normals were determined and the grids with similar surface normals were linked using a connected component analysis. Alharthy and Bethel [16] proposed a technique of moving windows to determine the surface slopes and then the planar surface segments were extracted. Tarsha-Kurdi et al. [17] applied Hough transform with RANSAC to determine plane surfaces of the roofs. Rabbani et al. [18] reviewed methods for modeling various surface types such as plane and curved surfaces (i.e., cylinder and sphere) from LiDAR data.Sohn and Dowman [19] used pan-sharpened multispectral images of IKONOS for building detection. The multispectral imagery provides useful information to identify buildings by utilizing land cover classification and image matching. However, there is limitation in use of imagery because of image quality such as low contrast, radiometric artifacts, occlusions due to central projection of the optical sensors, relief displacement, and shadows causing false breakline detection on the images. In addition, other limitations are the fact that both LiDAR data and imagery are not always available, and time difference between LiDAR data and image collection may result in undesirable outcomes due to possible surface change (e.g., construction of new buildings or earth work).Some methods utilize other geospatial sources (e.g., various imagery, maps, architectural drawings, and ground plans) to segment LiDAR data for 3D building modeling. Kim et al. [20] integrated stereo pair of the aerial images and LiDAR data for building detection and modeling. They derived 3D building boundaries by matching on the planar surface patches defined by the LiDAR data. Awrangjeb et al. [21] generated building models with masks derived from LiDAR data and multispectral images. Xiao et al. [22] implemented image matching algorithm for building modeling using LiDAR data and aerial images. Stilla and Jurkiewicz [23] used the building layer of the large-scale vector maps to separate and group parts of the buildings. They analyzed characteristic histograms of the laser altimeter data to identify roof type and to determine slope and orientation of the buildings for roof reconstruction. Vosselman and Dijkman [24] proposed a method for constructing 3D building models using 3D Hough transform and ground plans of the buildings to extract planar faces from the point clouds. The ground plans were utilized to segment planar faces of the roofs. Schenk and Csathό [25] proposed an approach to establish a common reference frame between LiDAR data and aerial images by utilizing sensor invariant features (i.e., edges and surface patches) for richer and more complete description of the surfaces. Thus, feature-level fusion of the LiDAR data and aerial imagery could contribute to mutual complement for each other’s shortcomings.Chang et al. [26] developed occlusion-based procedure for LiDAR data classification, building detection, and digital building model generation. The underlying concept of this method is based on a global operation to detect occluding points which belong to nonground objects. Wang et al. [27] presented automatic algorithm for building boundary detection. This approach includes identification of the building boundary by utilizing height and shape information for segmenting LiDAR data, boundary reconstruction using Hough transformation, and sequential linking. You and Lin [28] reconstructed building models and analyzed quality of the models by integrating LiDAR data and topographic maps. The limitations of integrating different sources of data are that all data should be available with temporal consistency. Cheng et al. [29] integrated multiview aerial imagery and LiDAR data to reconstruct 3D building models with accurate geometric position and high level of detail (LoD). The key tasks of the approach include determination of principal orientations of the buildings to improve boundary extraction, identification of boundary segments with K-means clustering, determination of the 3D boundary segments from LiDAR data and multiview images, and reconstruction of building models by recovering lost boundaries and rooftop patches using RANSAC algorithm. Rottensteiner and Clode [30] introduced a probabilistic context-based classification technique for building extraction and modeling. The procedure consists of semantic labelling of each LiDAR point using feature vectors, classification, and analysis of neighboring points for graph structures. Reconstruction of the building models was performed by estimation of model parameters and regularization.In recent years, deep learning (DL) based on the artificial neural networks is implemented to extract buildings and recognize 3D building models from LiDAR data [31–33]. Currently, most of the DL approaches in the field of remote sensing and geospatial engineering are focused on object detection, classification, and semantic segmentation primarily using imagery. In particular, majority of the works on object detection using DL apply convolutional neural networks. Despite the promising experiments performed, constructing optimal architecture of the neural network and exploring appropriate and effective learning algorithm for a specific purpose are main issues in DL [34–36]. It is expected that the DL approach will be extended to reconstruct various types of the 3D building models by utilizing variety of the geospatial data such as airborne and/or terrestrial LiDAR, multi/hyperspectral imagery, thermal-IR imagery, topographic maps, and various information derived from LiDAR data. Maltezos et al. [37] proposed multidimensional feature vector that consists of the raw LiDAR data, and additional features including entropy, height variation, intensity, distribution of normal vectors, number of returns, planarity, and standard deviation. The feature vector was used for training data with deep convolutional neural networks to extract buildings. The additional feature data was derived from the LiDAR data. Each feature data reflects unique physical property of the object.The segmented surface patches are used to determine model key features (e.g., building boundaries, intersecting lines or structure lines (e.g., ridges, hip), vertices, corners, apexes, and junctions) of the roofs. However, 3D building modeling with segmented surfaces requires topological relationship (i.e., adjacency) between surface patches. The proposed method aims to detect model key points (MKPs) of the buildings directly from the shaded relief images (SRIs) created from LiDAR data without segmentation. We utilized SRIs created from different locations of the light source (i.e., multidirectional SRIs) and then integrated MKPs detected from each SRI. For the better results, vertical exaggeration was performed to the LiDAR data before creating SRIs, and then anisotropic diffusion was carried out before detecting MKPs. The MKPs were detected by corner detector. Corner detection is frequently used in computer vision for motion detection, stereo image matching, image registration and mosaicking, object recognition, and modeling. The main objective of the corner detection is to determine distinctive, essential, and representative points of interest that could depict 3D geometric characteristics of the building such as shape and size. Eventually, automation of the digitizing (or feature collection) task for reconstruction of the building models can be achieved by detecting the corner points that are regarded as MKPs.There are considerable corner detectors available and some of well-known methods are Moravec, Förstner, Harris, smallest univalue segment assimilating nucleus (SUSAN), curvature scale space (CSS), scale invariant features transform (SIFT), speed-up robust features (SURF), features from accelerated segment test (FAST), and so on [38–46]. This paper implemented CSS algorithm to determine MKPs of the buildings with various roof types. The key features of the CSS are that, robust with respect to noise and superior to other detectors, the corner points are tracked through multiple scales to improve localization, false corners are removed, and T-junction corners can be detected. Detail description and performance of the CSS can be found in Mokhtarian and Suomela [42] and Zhang et al. [43]. In particular, capability of detecting T-junctions is important for reconstructing various roof types (e.g., hip, pyramid, and mansard). The proposed method was applied to various datasets. The results of the MKPs for each dataset were evaluated. The results from experiments demonstrate our method is quite efficient and practical for automatic determination of the MKPs of the buildings.


## 2. Proposed Methodology

In conventional photogrammetry, the MKPs are visually identified to measure coordinates of the MKPs by 3D digitizing on the stereo images. The main purpose of the proposed method is to determine locations of the MKPs of the buildings automatically using LiDAR data. The following procedures are commonly required in most of the methods for reconstructing 3D building models with LiDAR data:(i)
Separation of the point clouds that belong to buildings from other objects (e.g., ground surface, trees), so called filtering(ii)
Segmentation of the surface patches that form shape of the roofs based on the geometric similarity (e.g., slope and orientation of the slope)(iii)
Representation of each segmented roof surface with best-fitting mathematical function(iv)
Establishment of the adjacency relationship of the segmented surfaces to determine boundaries or intersections between surface patches.The concept behind the proposed method is shape from shading that is one of the computer vision problems for shape recognition and reconstruction. There are various visual cues to recognize objects. Shading is an important aspect of perceiving shape of the objects because different parts of the surface are oriented differently and thus appear with different brightness. The spatial variation of the brightness due to surface change (i.e., slope and orientation) is referred to as shading. Shading provides visual depth cue that makes possible appreciation of the surface topography and geometry [47]. In this regard, we utilized shaded relief images generated from LiDAR point clouds to detect more reliable and robust MKPs of the buildings.The essential procedures involved with reconstructing 3D building models with LiDAR data are building extraction and primitive (or feature, element) determination to depict building models. This paper proposes a method to determine MKPs that are one of the important elements to reconstruct building models. The corner detection algorithm used in digital image processing was applied to determine MKPs. Corner detection has been broadly used in motion detection, image registration, image matching, object recognition, and 3D modeling. The corners not only mean corner points themselves but also include vertices, apexes, and junctions formed by intersections of the edges or surfaces of the roof. Therefore, such features could be well-defined MKPs that are important points to define shapes of the individual buildings.The proposed method consists of the following procedures:(i)
Building extraction with marker-controlled watershed algorithm(ii)
Generation of multidirectional SRIs with vertically exaggerated data(iii)
Anisotropic diffusion to reduce noise effect of the SRIs(iv)
MKP determination using multidirectional SRIs by corner detector.


### 2.1. Building Extraction

Building extraction is an essential preprocessing because the goal is to determine MKPs of the buildings automatically. We applied marker-controlled watershed algorithm to extract building data. The marker-controlled watershed was developed to separate specific features from others by minimizing oversegmentation problem [48]. The concept of the watershed algorithm is to find catchment basins and watershed ridge lines by treating the LiDAR data as a terrain surface where light pixels are high and dark pixels are low. We improved the marker-controlled watershed method to avoid oversegmentation problem not only in the buildings but also on the ground surfaces. Figure 1 shows workflow of the building extraction based on the marker-controlled watershed.Figure 1
Workflow of building extraction scheme with marker-controlled watershed algorithm.A problem occurs due to imperfection of both foreground marker and background marker generation. This problem was resolved in the proposed method. It is crucial to determine building boundaries precisely in order for the marker-controlled watershed to work properly. If gradients by applying the Sobel operator are small (i.e., small slope), the buildings have flat roofs. Therefore, it is possible to extract buildings with flat roofs by setting threshold of the gradient (ε1). However, buildings have various roof types as well as flat roofs. The standard deviations of the moving window were computed to extract buildings with various roof types by setting threshold of the standard deviation (ε2) (see Algorithm 1).Algorithm 1:Improvement of mask in marker-controlled watershed.
Pgrad ← Modified slope by proposed method R, C ← Size of Mgrad (Modified slope) ε 1 , ε 2 ← Thresholds Mwin ← Moving window minV← Minimum value of Mgrad Pgrad = Mgrad Pgrad (Mgrad < ε1) = minV for y = 2 toR - 1 for x = 2 to C - 1 Mwin = Mgrad (n x nwindow at (x, y)) if (standard deviation of Mwin) < ε2 Pgrad (n x n window at (x, y)) = minV end end endFigure2 illustrates results from regular watershed and marker-controlled watershed. Marker-controlled watershed could provide quite useful results to detect MKPs from individually segmented building objects. Figure 8 shows building extraction results of the test datasets.Figure 2
Examples of watershed results. The data used to demonstrate watershed in this figure is part of Dongtan City. (a) LiDAR data. (b) Oversegmentation by regular watershed. (c) Segmentation by conventional marker-controlled watershed. (d) Segmentation by improved marker-controlled watershed.
(a) (b) (c) (d)


### 2.2. Generation of Shaded Relief Image

DSM and virtual light source are required to create SRIs. DSMs were generated using LiDAR data. Ground sampling distance (GSD) and accuracy of the DSM, location and direction of the light source, affect quality of the SRIs. The resolution of the SRI corresponds to the GSD of the DSM. Each pixel value of the SRIs is proportional to the amount of the light reflected from the surface element. Therefore, SRIs are determined by slope and aspect of the surface elements that are computed from DSM (see Figure3).Figure 3
Shaded relief image generation using DSM.The magnitude of the reflected light from each surface element to the pixel of the SRI is computed using(1)R=cos⁡I=cos⁡90°-βcos⁡S+sin⁡90°-βsin⁡Scos⁡α-A,whereR denotes magnitude of the reflected light to a pixel of the SRI,I denotes reflectance angle with respect to the surface normal of a surface element, α and β are azimuth and elevation angle of the light source, respectively, andS andA are slope and aspect of a surface element, respectively.Since SRI is imagery, the pixel value (i.e., brightness value) of the SRI is computed by(2)ISR=Gmax∙RwithGmax≤2n-1ifR<0thenR=0,where ISR is brightness of the SRI, Gmax denotes the maximum value ofn-bit image (e.g., 255 for 8-bit image), andn is number of the bit.The slope and aspect (i.e., orientation of the slope) can be derived from 3D coordinates of the 3x3 grid window of the DSM (see Figure4) using (4), (6), (7), and (8).Figure 4
3x3 grid window of DSM used for computing slope and aspect. X, Y, and Z are 3D coordinates of DSM.Slopes in X-direction are computed as follows:(3)∆X1=ZX+1,Y+1-ZX-1,Y+12∙GSDx∆X2=ZX+1,Y-ZX-1,Y2∙GSDx∆X3=ZX+1,Y-1-ZX-1,Y-12∙GSDx.Average slope in X-direction is(4)∆X=∆X1+∆X2+∆X33.Slopes in Y-direction are computed as follows:(5)∆Y1=ZX-1,Y+1-ZX-1,Y-12∙GSDy∆Y2=ZX,Y+1-ZX,Y-12∙GSDy∆Y3=ZX+1,Y+1-ZX+1,Y-12∙GSDy.Average slope in Y-direction is(6)∆Y=∆Y1+∆Y2+∆Y33.Therefore, slope and aspect at the center of the 3x3 grid window are (7)Slope≔tan-1∆X2+∆Y2(8)Aspect≔tan-1∆X∆Y.


### 2.3. Shaded Relief Image with Vertical Exaggeration

Figure5 demonstrates the effect of the vertical exaggeration to the SRI. The adjacent surfaces could be clearly distinguishable in the SRI by using vertically exaggerated DSM because vertical exaggeration changes surface slopes drastically.Figure 5
Shaded relief images before and after vertical exaggeration. (a) Before vertical exaggeration of LiDAR data. (b) After exaggeration of LiDAR data. (c) Shaded relief image before exaggeration. (d) Shaded relief image after exaggeration.
(a) (b) (c) (d)


### 2.4. Multidirectional Shaded Relief Images

As shown in Figure6, SRIs are differently generated depending on the location of the light source. In most cases, the NW location of the light sources is often used to improve visual perception for representing DSM or topographical maps. Our objective utilizing the multidirectional SRIs is to obtain more reliable MKPs by integrating results from each SRI. Integration of the multidirectional SRIs could provide 3D modeling with high LoD.Figure 6
Shaded relief images generated from different locations of the light source. (a) Shaded relief image with NW light source. (b) Shaded relief image with NE light source. (c) Shaded relief image with SW light source. (d) Shaded relief image with SE light source.


### 2.5. Anisotropic Diffusion

Reducing noise leads to better results in MKP detection. Anisotropic diffusion is a technique aiming at reducing or smoothing noises in the interior of the homogeneous regions without excessively removing or blurring significant parts of the object features such as edges, boundaries, or other details that are important for main processing. In this regard, anisotropic diffusion was implemented to effectively minimize local variation in the SRIs due to the vertical exaggeration.Edges and boundaries are preserved while noises are reduced to be homogeneous iteratively using (9), (10), and (11). More details of the anisotropic diffusion can be found in [49, 50].(9)∂Zx,y∂tt=∂Z2x,y∂x2t+∂Z2x,y∂y2t=div⁡cx,y∇Zx,yt(10)Zx,yt+1=Zx,y+T∑d=1n(cx,y∇Zdx,yt(11)cx,y=exp⁡∇Zx,yk2,whereZ(x, y) is LiDAR data,t denotes different level of the pyramid created by Gaussian filter, div stands for divergence operator,c(x, y) is diffusion function, ∇Z(x, y) is gradient,T is constant with respect to diffusion directions,n denotes number of directions (i.e., 4 or 8 directions at current center point), andk is diffusion coefficient.


### 2.6. Corner Detection and Determination of MKPs

Corners were detected using CSS algorithm. CSS is known as scale, rotation, and translation-invariant but sensitive to noises. CSS starts with edge detection with Canny operator. Then, gaps of the edges are connected (i.e., gap filling), and T-junctions are identified. Curvature is computed to determine corner candidates (i.e., maximum curvature provides provide corners). Figure7 shows CSS algorithm. The corner points are repeatedly traced through the scale space to improve positional accuracy of the corner points ([36, 37]).Figure 7
Workflow of corner detection with curvature scale space.Figure 8
Determination of a model key point by combining multidirectional shaded relief images. This figure illustrates determination of a model key point (i.e., apex) of the pyramid roof.MKPs detected from each SRI were combined by taking average coordinates of the MKPs. If MKP that is out of range, the specific point is excluded from the MKPs as shown in Figure8. The range value was six times GSD for each dataset in the experiments.


## 2.1. Building Extraction

Building extraction is an essential preprocessing because the goal is to determine MKPs of the buildings automatically. We applied marker-controlled watershed algorithm to extract building data. The marker-controlled watershed was developed to separate specific features from others by minimizing oversegmentation problem [48]. The concept of the watershed algorithm is to find catchment basins and watershed ridge lines by treating the LiDAR data as a terrain surface where light pixels are high and dark pixels are low. We improved the marker-controlled watershed method to avoid oversegmentation problem not only in the buildings but also on the ground surfaces. Figure 1 shows workflow of the building extraction based on the marker-controlled watershed.Figure 1
Workflow of building extraction scheme with marker-controlled watershed algorithm.A problem occurs due to imperfection of both foreground marker and background marker generation. This problem was resolved in the proposed method. It is crucial to determine building boundaries precisely in order for the marker-controlled watershed to work properly. If gradients by applying the Sobel operator are small (i.e., small slope), the buildings have flat roofs. Therefore, it is possible to extract buildings with flat roofs by setting threshold of the gradient (ε1). However, buildings have various roof types as well as flat roofs. The standard deviations of the moving window were computed to extract buildings with various roof types by setting threshold of the standard deviation (ε2) (see Algorithm 1).Algorithm 1:Improvement of mask in marker-controlled watershed.
Pgrad ← Modified slope by proposed method R, C ← Size of Mgrad (Modified slope) ε 1 , ε 2 ← Thresholds Mwin ← Moving window minV← Minimum value of Mgrad Pgrad = Mgrad Pgrad (Mgrad < ε1) = minV for y = 2 toR - 1 for x = 2 to C - 1 Mwin = Mgrad (n x nwindow at (x, y)) if (standard deviation of Mwin) < ε2 Pgrad (n x n window at (x, y)) = minV end end endFigure2 illustrates results from regular watershed and marker-controlled watershed. Marker-controlled watershed could provide quite useful results to detect MKPs from individually segmented building objects. Figure 8 shows building extraction results of the test datasets.Figure 2
Examples of watershed results. The data used to demonstrate watershed in this figure is part of Dongtan City. (a) LiDAR data. (b) Oversegmentation by regular watershed. (c) Segmentation by conventional marker-controlled watershed. (d) Segmentation by improved marker-controlled watershed.
(a) (b) (c) (d)


## 2.2. Generation of Shaded Relief Image

DSM and virtual light source are required to create SRIs. DSMs were generated using LiDAR data. Ground sampling distance (GSD) and accuracy of the DSM, location and direction of the light source, affect quality of the SRIs. The resolution of the SRI corresponds to the GSD of the DSM. Each pixel value of the SRIs is proportional to the amount of the light reflected from the surface element. Therefore, SRIs are determined by slope and aspect of the surface elements that are computed from DSM (see Figure3).Figure 3
Shaded relief image generation using DSM.The magnitude of the reflected light from each surface element to the pixel of the SRI is computed using(1)R=cos⁡I=cos⁡90°-βcos⁡S+sin⁡90°-βsin⁡Scos⁡α-A,whereR denotes magnitude of the reflected light to a pixel of the SRI,I denotes reflectance angle with respect to the surface normal of a surface element, α and β are azimuth and elevation angle of the light source, respectively, andS andA are slope and aspect of a surface element, respectively.Since SRI is imagery, the pixel value (i.e., brightness value) of the SRI is computed by(2)ISR=Gmax∙RwithGmax≤2n-1ifR<0thenR=0,where ISR is brightness of the SRI, Gmax denotes the maximum value ofn-bit image (e.g., 255 for 8-bit image), andn is number of the bit.The slope and aspect (i.e., orientation of the slope) can be derived from 3D coordinates of the 3x3 grid window of the DSM (see Figure4) using (4), (6), (7), and (8).Figure 4
3x3 grid window of DSM used for computing slope and aspect. X, Y, and Z are 3D coordinates of DSM.Slopes in X-direction are computed as follows:(3)∆X1=ZX+1,Y+1-ZX-1,Y+12∙GSDx∆X2=ZX+1,Y-ZX-1,Y2∙GSDx∆X3=ZX+1,Y-1-ZX-1,Y-12∙GSDx.Average slope in X-direction is(4)∆X=∆X1+∆X2+∆X33.Slopes in Y-direction are computed as follows:(5)∆Y1=ZX-1,Y+1-ZX-1,Y-12∙GSDy∆Y2=ZX,Y+1-ZX,Y-12∙GSDy∆Y3=ZX+1,Y+1-ZX+1,Y-12∙GSDy.Average slope in Y-direction is(6)∆Y=∆Y1+∆Y2+∆Y33.Therefore, slope and aspect at the center of the 3x3 grid window are (7)Slope≔tan-1∆X2+∆Y2(8)Aspect≔tan-1∆X∆Y.


## 2.3. Shaded Relief Image with Vertical Exaggeration

Figure5 demonstrates the effect of the vertical exaggeration to the SRI. The adjacent surfaces could be clearly distinguishable in the SRI by using vertically exaggerated DSM because vertical exaggeration changes surface slopes drastically.Figure 5
Shaded relief images before and after vertical exaggeration. (a) Before vertical exaggeration of LiDAR data. (b) After exaggeration of LiDAR data. (c) Shaded relief image before exaggeration. (d) Shaded relief image after exaggeration.
(a) (b) (c) (d)


## 2.4. Multidirectional Shaded Relief Images

As shown in Figure6, SRIs are differently generated depending on the location of the light source. In most cases, the NW location of the light sources is often used to improve visual perception for representing DSM or topographical maps. Our objective utilizing the multidirectional SRIs is to obtain more reliable MKPs by integrating results from each SRI. Integration of the multidirectional SRIs could provide 3D modeling with high LoD.Figure 6
Shaded relief images generated from different locations of the light source. (a) Shaded relief image with NW light source. (b) Shaded relief image with NE light source. (c) Shaded relief image with SW light source. (d) Shaded relief image with SE light source.


## 2.5. Anisotropic Diffusion

Reducing noise leads to better results in MKP detection. Anisotropic diffusion is a technique aiming at reducing or smoothing noises in the interior of the homogeneous regions without excessively removing or blurring significant parts of the object features such as edges, boundaries, or other details that are important for main processing. In this regard, anisotropic diffusion was implemented to effectively minimize local variation in the SRIs due to the vertical exaggeration.Edges and boundaries are preserved while noises are reduced to be homogeneous iteratively using (9), (10), and (11). More details of the anisotropic diffusion can be found in [49, 50].(9)∂Zx,y∂tt=∂Z2x,y∂x2t+∂Z2x,y∂y2t=div⁡cx,y∇Zx,yt(10)Zx,yt+1=Zx,y+T∑d=1n(cx,y∇Zdx,yt(11)cx,y=exp⁡∇Zx,yk2,whereZ(x, y) is LiDAR data,t denotes different level of the pyramid created by Gaussian filter, div stands for divergence operator,c(x, y) is diffusion function, ∇Z(x, y) is gradient,T is constant with respect to diffusion directions,n denotes number of directions (i.e., 4 or 8 directions at current center point), andk is diffusion coefficient.


## 2.6. Corner Detection and Determination of MKPs

Corners were detected using CSS algorithm. CSS is known as scale, rotation, and translation-invariant but sensitive to noises. CSS starts with edge detection with Canny operator. Then, gaps of the edges are connected (i.e., gap filling), and T-junctions are identified. Curvature is computed to determine corner candidates (i.e., maximum curvature provides provide corners). Figure7 shows CSS algorithm. The corner points are repeatedly traced through the scale space to improve positional accuracy of the corner points ([36, 37]).Figure 7
Workflow of corner detection with curvature scale space.Figure 8
Determination of a model key point by combining multidirectional shaded relief images. This figure illustrates determination of a model key point (i.e., apex) of the pyramid roof.MKPs detected from each SRI were combined by taking average coordinates of the MKPs. If MKP that is out of range, the specific point is excluded from the MKPs as shown in Figure8. The range value was six times GSD for each dataset in the experiments.


## 3. Experimental Results and Discussion

In this section, description of the test datasets, experimental results, and analysis of the results are presented.


### 3.1. Test Datasets

We evaluate the proposed method on five airborne LiDAR datasets including simulated and real datasets (see Figure9). The simulated data consists of various roof types (e.g., gable, hip, half-hip, pyramid, gambrel, saw-tooth, and dome) with 0.25 m GSD. The real data includes urban area (Dongtan City) in Korea (0.50 m GSD), part of the University of Calgary (0.60 m GSD), and Vaihingen data (0.50 m GSD). Vaihingen data is a benchmark dataset provided by the German Society for Photogrammetry, Remote Sensing, and Geoinformation (DGPF) [44]. All datasets were resampled into regular grids using the nearest neighbor interpolation.Figure 9
Experimental test datasets. (a) Simulated building models. (b) Simulated LiDAR data. (c) Aerial photo of Dongtan City. (d) Dongtan City LiDAR data. (e) Aerial photo of University of Calgary. (f) University of Calgary LiDAR data. (g) Color-IR orthoimage of Vaihingen data. (h) Vaihingen LiDAR data.
(a) (b) (c) (d) (e) (f) (g) (h)


### 3.2. Building Extraction

Table1 shows the parameter setting for experiments of the building extraction. Figure 10 represents results of the building extraction for each dataset. The results show that most of the buildings were properly extracted based on our observation.Table 1
Parameters and corresponding thresholds for building extraction.
Procedure Parameter Threshold value Marker-controlled watershed FBM Morphological filter[pixel] 5 × 5 BGM Binarization[m] Mode + 5 Modified mask Slope (ε1) [%] 200 Std. dev (ε2) [%] ±20  Extracting Building data Height[m] Mode + 3 Area[m2] 75 < Area < 3,750 Std. dev[m] ±1Figure 10
Extracted buildings by marker-controlled watershed algorithm. (a) Simulated data. (b) Dongtan City data. (c) University of Calgary data. (d) Vaihingen data.
(a) (b) (c) (d)


### 3.3. Corner Detection and MKPs

Different corner detectors (Moravec, Harris, FAST, and CSS) were evaluated using simulated dataset. The CSS is superior to other corner detectors based on our experiments and observation as shown in Figure11. Most corner points were detected by CSS while some important points were not detected by other methods. In this regard, CSS was applied to determine MKPs in this study. It is noticed that raw LiDAR data was used to evaluate corner detectors. In other words, SRIs were not used for this evaluation test but the raw LiDAR point clouds with 3D coordinates (i.e., LiDAR depth image) were used.Figure 11
Performance of the corner detectors. (a) Moravec. (b) Harris. (c) FAST. (d) CSS.
(a) (b) (c) (d)Table2 shows the parameters applied to the experiments for MKP determination. Results of the determined MKPs for each dataset are represented from Figures 12–18. The results include MKPs from multidirectional SRIs (see Figures 12, 14, 16, and 18), and a comparison between MKPs determined by proposed approach and from raw data; Figures 13, 15, 17, and 19 show a comparison between results from the proposed method (see Figure (a)s), and from raw LiDAR data without any further processing (see Figure (b)s).Table 2
Parameters and corresponding values for experiment to determine MKPs.
Procedure Parameter Value Location of light source Azimuth[​°] 45 (NE), 135 (SE), 225 (SW), 315 (NW) Elevation[​°] 45  Vertical exaggeration Z-scale factor 3  CSS Canny threshold Low: 0.1, High: 0.2 Curvature 0.2 σ 3 Gap distance [pixel] 2  Anisotropic diffusion Number of iterations 20 T 0.2 k 3Figure 12
MKPs from each SRI of simulated data. (a) NW direction. (b) NE direction. (c) SW direction. (d) SE direction.
(a) (b) (c) (d)Figure 13
MKPs of simulated data. (a) Result from the proposed method. (b) Result from raw LiDAR data.
(a) (b)Figure 14
MKPs from each SRI of Dongtan City data. (a) NW direction. (b) NE direction. (c) SW direction. (d) SE direction.
(a) (b) (c) (d)Figure 15
MKPs of Dongtan City data. (a) Result from the proposed method. (b) Result from raw LiDAR data.
(a) (b)Figure 16
MKPs from each SRI of University of Calgary data. (a) NW direction. (b) NE direction. (c) SW direction. (d) SE direction.
(a) (b) (c) (d)Figure 17
MKPs of University of Calgary data. (a) Result from the proposed method. (b) Result from raw LiDAR data.
(a) (b)Figure 18
MKPs from each SRI of Vaihingen data. (a) NW direction. (b) NE direction. (c) SW direction. (d) SE direction.
(a) (b) (c) (d)Figure 19
MKPs of Vaihingen data. (a) Result from the proposed method. (b) Result from raw LiDAR data.
(a) (b)Figures13, 15, 17, and 19(a) show the final MKPs of the buildings after combining multidirectional SRI with vertical exaggeration and anisotropic diffusion. On the other hand, Figures 13, 15, 17, and 19(b) show the MKPs by applying CSS corner detector to the raw LiDAR data without any further processing. As shown in the figures, most of the important MKPs could be determined by applying the proposed method including apexes of the pyramid roofs. However, applying the corner detector to the raw data results in detecting corners of the buildings. Therefore, the proposed method would be effective and beneficial to construct 3D building models with high LoD.


### 3.4. Analysis and Discussion

The results from each procedure involved with proposed approach were analyzed. To evaluate performance of our approach, we assessed both detection and undetection rates, and false detection rates for each dataset. The results show that improvement has been achieved by the proposed method. Overall quality of the results was also provided (see Tables3 and 4). The total number of true MKPs was identified by visual inspection on the aerial imagery. It is noticed that the “false detected MKP” in Tables 3 and 4 is not actually MPKs, but points detected by the corner detector. For the comparison purpose, Tables 3 and 4 show evaluation results using raw LiDAR data and results from the proposed method, respectively.Table 3
Evaluation of MKP detection from raw LiDAR data.
Dataset Total number of true MKPs Detected MKP Undetected MKP Falsely detected MKP No. Rate[%] No. Rate[%] No. Rate[%] Simulation 150 82 55 68 45 0 0 Dongtan City 203 108 53 95 47 6 3 University of Calgary 132 88 67 44 33 6 5 Vaihingen 218 138 63 80 37 0 0  Overall 703 416 59 287 41 12 2Table 4
Evaluation of MKP detection from proposed method.
Dataset Total number of true MKPs Detected MKP Undetected MKP Falsely detected MKP No. Rate[%] No. Rate[%] No. Rate[%] Simulation 150 118 79 32 21 19 13 Dongtan City 203 140 69 63 31 27 13 University of Calgary 132 98 74 34 26 39 30 Vaihingen 218 140 64 78 36 4 2  Overall 703 496 71 207 29 89 13Comparing the results from our method with results from raw LiDAR data, overall accuracy increases by 12% in terms of detection rate. However, false alarm rate (i.e., falsely detected MKPs) increases by 11%. This may be caused from vertical exaggeration and interpolation of the LiDAR data even though anisotropic diffusion to reduce noises was applied. In case of the Vaihingen data, there is no improvement in terms of detection rate (i.e., number of points detected). However, the MKPs required for high LoD modeling were detected by the proposed method. When the raw LiDAR data was used, the MKPs were not properly detected. Point density (i.e., GSD) and interpolation of the LiDAR point clouds might influence the accuracy of the results. In particular, the resolution (i.e., GSD) of Vaihingen orthoimage is 0.09 m while GSD of the LiDAR data is 0.50 m. Therefore, the detection rate of our method might be underestimated due to the extreme resolution difference between ground truth and LiDAR data. Majority of the false detection of MKPs are found along building boundaries due to interpolation of the LiDAR data. In consequence, refinement of the building boundaries by regularization with building hypotheses and/or constraints could solve false detection of the MKPs.


## 3.1. Test Datasets

We evaluate the proposed method on five airborne LiDAR datasets including simulated and real datasets (see Figure9). The simulated data consists of various roof types (e.g., gable, hip, half-hip, pyramid, gambrel, saw-tooth, and dome) with 0.25 m GSD. The real data includes urban area (Dongtan City) in Korea (0.50 m GSD), part of the University of Calgary (0.60 m GSD), and Vaihingen data (0.50 m GSD). Vaihingen data is a benchmark dataset provided by the German Society for Photogrammetry, Remote Sensing, and Geoinformation (DGPF) [44]. All datasets were resampled into regular grids using the nearest neighbor interpolation.Figure 9
Experimental test datasets. (a) Simulated building models. (b) Simulated LiDAR data. (c) Aerial photo of Dongtan City. (d) Dongtan City LiDAR data. (e) Aerial photo of University of Calgary. (f) University of Calgary LiDAR data. (g) Color-IR orthoimage of Vaihingen data. (h) Vaihingen LiDAR data.
(a) (b) (c) (d) (e) (f) (g) (h)


## 3.2. Building Extraction

Table1 shows the parameter setting for experiments of the building extraction. Figure 10 represents results of the building extraction for each dataset. The results show that most of the buildings were properly extracted based on our observation.Table 1
Parameters and corresponding thresholds for building extraction.
Procedure Parameter Threshold value Marker-controlled watershed FBM Morphological filter[pixel] 5 × 5 BGM Binarization[m] Mode + 5 Modified mask Slope (ε1) [%] 200 Std. dev (ε2) [%] ±20  Extracting Building data Height[m] Mode + 3 Area[m2] 75 < Area < 3,750 Std. dev[m] ±1Figure 10
Extracted buildings by marker-controlled watershed algorithm. (a) Simulated data. (b) Dongtan City data. (c) University of Calgary data. (d) Vaihingen data.
(a) (b) (c) (d)


## 3.3. Corner Detection and MKPs

Different corner detectors (Moravec, Harris, FAST, and CSS) were evaluated using simulated dataset. The CSS is superior to other corner detectors based on our experiments and observation as shown in Figure11. Most corner points were detected by CSS while some important points were not detected by other methods. In this regard, CSS was applied to determine MKPs in this study. It is noticed that raw LiDAR data was used to evaluate corner detectors. In other words, SRIs were not used for this evaluation test but the raw LiDAR point clouds with 3D coordinates (i.e., LiDAR depth image) were used.Figure 11
Performance of the corner detectors. (a) Moravec. (b) Harris. (c) FAST. (d) CSS.
(a) (b) (c) (d)Table2 shows the parameters applied to the experiments for MKP determination. Results of the determined MKPs for each dataset are represented from Figures 12–18. The results include MKPs from multidirectional SRIs (see Figures 12, 14, 16, and 18), and a comparison between MKPs determined by proposed approach and from raw data; Figures 13, 15, 17, and 19 show a comparison between results from the proposed method (see Figure (a)s), and from raw LiDAR data without any further processing (see Figure (b)s).Table 2
Parameters and corresponding values for experiment to determine MKPs.
Procedure Parameter Value Location of light source Azimuth[​°] 45 (NE), 135 (SE), 225 (SW), 315 (NW) Elevation[​°] 45  Vertical exaggeration Z-scale factor 3  CSS Canny threshold Low: 0.1, High: 0.2 Curvature 0.2 σ 3 Gap distance [pixel] 2  Anisotropic diffusion Number of iterations 20 T 0.2 k 3Figure 12
MKPs from each SRI of simulated data. (a) NW direction. (b) NE direction. (c) SW direction. (d) SE direction.
(a) (b) (c) (d)Figure 13
MKPs of simulated data. (a) Result from the proposed method. (b) Result from raw LiDAR data.
(a) (b)Figure 14
MKPs from each SRI of Dongtan City data. (a) NW direction. (b) NE direction. (c) SW direction. (d) SE direction.
(a) (b) (c) (d)Figure 15
MKPs of Dongtan City data. (a) Result from the proposed method. (b) Result from raw LiDAR data.
(a) (b)Figure 16
MKPs from each SRI of University of Calgary data. (a) NW direction. (b) NE direction. (c) SW direction. (d) SE direction.
(a) (b) (c) (d)Figure 17
MKPs of University of Calgary data. (a) Result from the proposed method. (b) Result from raw LiDAR data.
(a) (b)Figure 18
MKPs from each SRI of Vaihingen data. (a) NW direction. (b) NE direction. (c) SW direction. (d) SE direction.
(a) (b) (c) (d)Figure 19
MKPs of Vaihingen data. (a) Result from the proposed method. (b) Result from raw LiDAR data.
(a) (b)Figures13, 15, 17, and 19(a) show the final MKPs of the buildings after combining multidirectional SRI with vertical exaggeration and anisotropic diffusion. On the other hand, Figures 13, 15, 17, and 19(b) show the MKPs by applying CSS corner detector to the raw LiDAR data without any further processing. As shown in the figures, most of the important MKPs could be determined by applying the proposed method including apexes of the pyramid roofs. However, applying the corner detector to the raw data results in detecting corners of the buildings. Therefore, the proposed method would be effective and beneficial to construct 3D building models with high LoD.


## 3.4. Analysis and Discussion

The results from each procedure involved with proposed approach were analyzed. To evaluate performance of our approach, we assessed both detection and undetection rates, and false detection rates for each dataset. The results show that improvement has been achieved by the proposed method. Overall quality of the results was also provided (see Tables3 and 4). The total number of true MKPs was identified by visual inspection on the aerial imagery. It is noticed that the “false detected MKP” in Tables 3 and 4 is not actually MPKs, but points detected by the corner detector. For the comparison purpose, Tables 3 and 4 show evaluation results using raw LiDAR data and results from the proposed method, respectively.Table 3
Evaluation of MKP detection from raw LiDAR data.
Dataset Total number of true MKPs Detected MKP Undetected MKP Falsely detected MKP No. Rate[%] No. Rate[%] No. Rate[%] Simulation 150 82 55 68 45 0 0 Dongtan City 203 108 53 95 47 6 3 University of Calgary 132 88 67 44 33 6 5 Vaihingen 218 138 63 80 37 0 0  Overall 703 416 59 287 41 12 2Table 4
Evaluation of MKP detection from proposed method.
Dataset Total number of true MKPs Detected MKP Undetected MKP Falsely detected MKP No. Rate[%] No. Rate[%] No. Rate[%] Simulation 150 118 79 32 21 19 13 Dongtan City 203 140 69 63 31 27 13 University of Calgary 132 98 74 34 26 39 30 Vaihingen 218 140 64 78 36 4 2  Overall 703 496 71 207 29 89 13Comparing the results from our method with results from raw LiDAR data, overall accuracy increases by 12% in terms of detection rate. However, false alarm rate (i.e., falsely detected MKPs) increases by 11%. This may be caused from vertical exaggeration and interpolation of the LiDAR data even though anisotropic diffusion to reduce noises was applied. In case of the Vaihingen data, there is no improvement in terms of detection rate (i.e., number of points detected). However, the MKPs required for high LoD modeling were detected by the proposed method. When the raw LiDAR data was used, the MKPs were not properly detected. Point density (i.e., GSD) and interpolation of the LiDAR point clouds might influence the accuracy of the results. In particular, the resolution (i.e., GSD) of Vaihingen orthoimage is 0.09 m while GSD of the LiDAR data is 0.50 m. Therefore, the detection rate of our method might be underestimated due to the extreme resolution difference between ground truth and LiDAR data. Majority of the false detection of MKPs are found along building boundaries due to interpolation of the LiDAR data. In consequence, refinement of the building boundaries by regularization with building hypotheses and/or constraints could solve false detection of the MKPs.


## 4. Conclusions

In this paper, we proposed a novel method to determine MKPs of the buildings using SRIs generated from airborne LiDAR data. The procedure includes building extraction, vertical exaggeration, anisotropic diffusion, corner detection, and MKP determination. The MKPs determined from the proposed approach provide efficient way to feasible results to reconstruct 3D building models. The quality of the results was analyzed and evaluated by comparing true MKPs with MKPs determined by our approach. The following aspects have been highlighted:(i)
Extracting building objects is an important and essential process for 3D building modeling. We have been able to achieve precise results by improving the existing marker-controlled watershed algorithm.(ii)
The corners detected in each SRI that is created from different location of the light source could mutually compensate for more complete detection of the MKPs.(iii)
Detail MKPs such as apexes of the pyramid roofs, intersections of the gable and hip roofs, and superstructures of the roofs can be detected.(iv)
We achieved an overall improvement of around 12%. However, falsely detected MKPs have increased. This problem is probably caused by noise emphasis due to height exaggeration. Since we expected this problem, the anisotropic diffusion was implemented, but the noises might not be completely removed.A future work of this study will be extended to reconstruct complete building modeling such as 3D wireframe building models by connecting the MKPs. This task might require topological analysis of the MKPs for appropriate connection of the MKPs that represent individual buildings.



---

*Source: 2985014-2019-04-21.xml*